











* ante 


on aa 





QUARTERLY OF APPLIED MATHEMATICS 


Vol. XVIII JANUARY, 1961 No. 4 








THE MISES YIELD CONDITION FOR 
ROTATIONALLY SYMMETRIC SHELLS* 


BY 
P. G. HODGE, Jr. 
Illinois Institute of Technology, Chicago, Illinois 


Summary. ‘The stress state in a rotationally symmetric shell under small displace- 
ments is characterized by the direct stresses and moments in the circumferential and 
longitudinal directions. If the shell material is perfectly plastic, it is desirable to express 
the material yield condition in terms of these stress resultants. Previous investigations 
have obtained this yield condition in certain special cases and for the maximum shear 
stress criterion in the general case. 

Here, a derivation is given based on the Mises or octahedral-shear-stress criterion 
for both uniform and idealized-sandwich shells. The flow law relating extension and 
curvature rates of the middle surface to the stress state is also obtained. The general 
equations are obtained in closed parametric form and various special cases are explicitly 
presented. 

1. Introduction. In the analysis of structures which are thin in one direction it is 
generally convenient to deal with stress resultants integrated over the thickness, rather 
than with the stresses themselves. The present paper will be directly concerned with 
rotationally symmetric shells. However, it will be shown that various simpler structures 
such as plates, slabs, arches, and frames may be regarded as suitable special cases of 
a theory based on rotationally symmetric shells. Although the methods used can theo- 
retically be applied to still more general shell problems, the difficulties are formidable 
and the subject will not be treated here. 

The state of stress in a symmetrically loaded rotationally symmetric shell is specified 
by circumferential and longitudinal direct stresses and moments and by a shear force. 
However, the basic assumptions are made that straight lines normal to the median 
surface of the shell remain straight and normal to the deformed median surface and 
that the displacements are small. It follows that shear strains are neglected so that 
the shear force is a reaction, not a generalized stress. Therefore, there are four generalized 
stresses, VN, , VN, , MW, ,M/, . These quantities are related to the physical stress com- 


ponents by 
aH 


H 
N.=[ o.dZ, M, - | Ze. dZ, (1.1) 


a 
-H “-H 


where the subscript a may stand for either @ or ¢. 


*Received August 10, 1959. This investigation was supported by the United States Office of Naval 
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As first shown by Prager [1] the generalized strain rates corresponding to these 
generalized stresses are the extension rates A, , 4, and curvature rates K, , K, of the 
middle surface of the shell. In terms of these quantities, the strain rate components 
at a point are 

é.=rA,+ ZK, . (1.2) 


In order to formulate the shell problem, it is necessary to express all of the constituent 
equations in terms of the generalized stresses and strain rates. The equations of equili- 
brium and the relations between strain rates and velocities are easily derived by direct 
consideration of an infinitesimal element of the shell (see, for example [2]). In the case 
of an elastic material, the remaining equations necessary for a solution are obtained 
by combining Hooke’s law with Eqs. (1.1) and (1.2) to obtain the elastic relations 
between generalized strain rate and stress. 

For a perfectly plastic material, it is first necessary to express the yield condition 
in terms of generalized stresses. For various particular cases, such as beams under 
combined tension and bending [3] or circular plates [4], this is easily done for any yield 
condition. However, for the general rotationally symmetric shell, the problem becomes 
more difficult. Yield conditions have been obtained for shells where the material satisfies 
Tresca’s yield condition of maximum shear [4-8]. In the present paper, we shall derive 
the yield condition for a shell whose material satisfies Mises yield condition of maximum 
octahedral stress. 

Once the yield condition has been derived, the remaining equations are furnished 
by the plastic potential flow law. In geometrical terms, this states that the strain rate 
vector whose components are the generalized strain rates must be normal to the yield 
surface. 

The present paper is concerned solely with the derivation of the yield condition 
and flow law for a rotationally symmetric shell. Applications of the theory will be 
reported on elsewhere [9]. 

2. Derivation of yield surface. As generalized stresses and strain rates for the 
rotationally symmetric shell we choose the dimensionless quantities 


ne =WN,/N,, m, =M./M,, Ries ka = (M,/No)K. (2.1) 


a 


where 


No = 2Ha ,; M, = H’s, (2.2) 


a» being the yield stress and 2H the shell thickness. The dimensionless counterparts 
of Eqs. (1.1) and (1.2) are then 
1 1 
n, = (1/2) [ 8, dz, mm. = [ 28, dz r 
Jay Ji (2.3) 
Eo = Ag + 2Wka ; z= Z/H, Sa = 02/00; 
At a plastic point the stresses must satisfy the Mises yield condition 

si — 8,8 +8; = 1 (2.4) 
and the associated plastic potential flow rule. This latter states that the strain rate 
vector (€, , €2) is normal to the curve (2.4), hence 


€, = v(2s, — 8), €2 = v(28, — 8), (2.5) 
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where v is an arbitrary positive scalar which describes the indeterminate magnitude 
of the strain rate vector. 

By combining Eqs. (2.3)-(2.5) we can establish the generalized stresses as homo- 
geneous functions of order zero of the generalized strain rates. It follows that the gener- 
alized stresses are thus dependent only upon three parameters representing the direction 
of the generalized strain rate vector in a four dimensional space. Since the representation 
of four variables in terms of three parameters is equivalent to a single equation relating 
these variables, we can thus obtain the desired yield condition. 

There are, of course, many ways of choosing three parameters to represent the 
direction of the strain-rate vector; the following will prove convenient for later inte- 
gration 

(A, , Ax) = +2rt cos (r + y + 2/6), 
(2.6) 
(x, , Ko) = +v cos (r + 7/6). 


Here, and in the following, the first and second terms in parenthesis are to be associated 
with the upper and lower signs, respectively. 
In terms of r, ¢t, and y, the stresses are 
2 tsin(r + y + 2/6) + zsin (r F 2/6) 
, 82) ee’; a * “yee \2 — ————— (2.7) 
[(z + t cos y) + (tsin y)"] 





(s 

l s" 
Substitution of Eqs. (2.7) into Eqs. (2.3) and integration will give the desired yield 
condition. To carry out this integration we make the substitution defined in Fig. 1 


and introduce new parameters p and q as the values of w at z = —1l andz = +1, re- 
spectively. If ¢ is finite and ¢ sin y # 0, then p and q must satisfy one or the other of 


—r/2<pswiiq<7/2, (2.8a) 
r/2<psw<q < 3n/2. (2.8b) 


In either case, the generalized stresses and strain rates are obtained in terms of p, q, 
and r [10]: 
boa (1 + sin g)(1 — sin p) 
(n; , n)[3'’* sin (p — g)] = 208 s (r / nT age Y : 
N, , N2)[3'’ sin (p — q)] = cos p cos q cos (r ¥ 1/6) log () —sin Xl + op 


+ 2sin (r + 1/6)(cos p — cos q) (2.9) 





(m, , m2)[3'”’ sin’ (p — q)/2 cos p cos qg)] 
= [sin (p + q) cos (r = +/6) + cos p cos gsin (r F 7/6)] 
-log [((1 + sin g)(1 — sin p)(1 — sin g) "(1 + sin p)'] 
— 4 cos (r * 2/6)(cos p — cos g) — 2sin (r F 7/6)(sin g — sin p) 
(A; , A2) sin (p — gq) = »o[42 cos gsin (r — p + 1/6) + 2 cos psin (r — gq + 7/6)] 
(xk, , Ko) = vy cos (r + 2/6) sin (p — Qq). 


Equations (2.9) give two different hypersurfaces, depending upon whether p and g 
satisfy Eq. (2.8a) or (2.8b). The dividing “hypercurve” between them may be found 
directly by setting sin y = 0 before integrating the generalized stresses [10] or by letting 
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z+! cosy 








tsin y 


Fia. 1. Definition of w 


p and/or g tend to +7/2 in Eqs. (2.9). The result is 
(n, , M2) = (2/3'”’)tsin (r F 7/6) 
(m, , M2) = (2/3'*)(1 — #) sin (r + 7/6) (2.10) 
-l1<i<l. 

Much simpler expressions are obtained for an ideal sandwich shell composed of 
two thin sheets of thickness J each, separated by a core of thickness 2H’. The sheets 
are so thin that the stress variation across each sheet can be neglected; the core has 
no tensile strength but can carry the necessary shear. Evidently Eqs. (2.2) and the 
last Eq. (2.3) should be replaced by 

No = 2ei/, M, = 20;H’' J, = ee. (2.11) 

It follows from statics that the stress resultants are related to the stresses s{ and s, 

in the top and bottom sheets, respectively, by 


So = Na tM, . = 4. ~— WM, . (2:32) 


a 


The yield condition (2.4) cannot be violated by the stresses in either sheet. Therefore, 
substitution of Eqs. (2.12) in (2.4) shows that the yield condition consists of the two 
non-linear surfaces 


(n, & m,)? — (n, & m,)(m2 + me) + (mn.  m.) = 1. (2.13a, b) 


The corresponding strain rates are 


A; = v [2(n, + m,) — (nm. + m,)], A> = v [2(nm2 + m.) — (n, + m)] 2 14a. b) 
(2. a, 0 
ky = ty [2(n, + m,) — (m2 + m.)], «2 = 4v°[2(n. + m) — (n, + m,)]. 


If the stresses satisfy both of Eqs. (2.13), then the strain rates may be any combination 
of those in Eqs. (2.14a and b), provided only that the coefficients v* and »” are both 
non-negative. 

3. Special cases. Various structures of practical importance may be regarded as 
special cases of shells of revolution. Included in this category are arches, circular slabs 
under rotationally symmetric in-plane loads, rotational bending of a circular plate, 
and rotationally loaded circular cylindrical shells. In each case, the yield condition 
can be obtained as an appropriate special case of Eqs. (2.9) or (2.13). The following 
results of interest are then obtained. 

a. Circular cylindrical shell. By assumption, xe = 0 hence m, = m, is a reaction 
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Fig. 2. Yield curves for circular cylindrical shell without end load. 
— —  —WMises condition, uniform shell 
cette eeee Mises condition, sandwich shell 
——-——Tresca condition, uniform shell 
—~ — ——Tresca condition, sandwich shell 


to be eliminated from the yield condition. For the sandwich shell this process leads 
to two surfaces: 
22n, — - = “ 271/2 ¥. 271/2 (3.1) 
+2(2n, — n,) = [4 — 3(n, m,) |" + [4 — 3(n, + m,)]. 3. 
For the uniform shell we obtain the two surfaces 
2 cos p — cos q 
il 


3°" sin(p — q) 


7 __ COS p COS q_ (1 + sin q) (1 — sin p) cos p — cos g | 
= | 2 sin (p — q) log (1 — sin g) (1 + sin p) + 3'” sin (p — gj’ 


n, + 


(3.2) 


m, = 


| 2 cos BP cos’ q lox (1 + sin 9) — sin p) 
3°" sin’ (p — q) (1 — sin g)(1 + sin p) 


, £08 p cos q (sin g — sing) |. 
3'” sin’ (p — q) 


b. Circular cylindrical shell without end load. This case is a special case of the 
preceding one obtained by setting n, = 0 in Eqs. (3.1) or (3.2). We obtain 


ns + (3/4)m2 = 1 (3.3) 
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Fig. 3. Yield curves for stresses or circular plate under pure bending or pure tension 
Mises condition, both shells 
Tresca condition, both shells 
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Fic. 4. Yield curves for arches 
Uniform shell, both conditions 
— - —Sandwich shell, both conditions 
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for the sandwich shell. For the uniform shell 


+(1/2) cot q log itees 


“ sitet (3.4) 
mm. = pa E t ql 1 + sin ¢ — 2 ese | 
™ wee rs | —sing ~ q 


c. Circular slab under tension. Here m, = mz = 0 by assumption. The entire yield 
curve for the uniform shell is on the locus (2.10) of singular points. Setting ¢ = 1 and 
eliminating r from the first line of (2.10) we obtain 


n> — nn. + nz = 1. (3.5) 
It follows from Eqs. (2.13) that (3.5) is also valid for the sandwich shell. 
d. Circular plate under bending. In this case n, = n, = 0. Reasoning similar to 
that in case (c) shows that the yield curve for the uniform or sandwich plate is 
m; — mm, + m; = 1. (3.6) 


e. Curved beam or arch. Here transverse stresses are negligible, so that n, = m, = 0. 
Here also the singular curve for the uniform shell applies in the form 


my = +(1 — n}), (3.7) 
whereas the sandwich shell reduces to the linear expressions 
Ne — mM = +1, Ne +m, = +1. (3.8) 


Figures 2, 3, and 4 show the special yield curves obtained in this section. For com- 


parison, the corresponding curves for materials which satisfy the Tresca condition are 


also shown. 
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BOOK REVIEWS 


On numerical approximation. Edited by Rudolph E. Langer. The University of Wisconsin 
Press, Madison, 1959. x + 462 pp. $4.50. 


This book presents the 21 papers which were delivered at the Symposium on Numerical Approxi- 
mation conducted by the Mathematics Research Center, United States Army, at the University of 
Wisconsin from April 21 to 23, 1958. The objective of this symposium was the presentation and dis- 
cussion of recent developments in the field of numerical approximation centered around three general 
themes: linear approximation, extremal approximation, and algorithms. 

In “On Trends and Problems in Numerical Approximations” A. M. Ostrowski surveys, enter- 
tainingly, developments in interpolation theory, in methods of successive approximation and in the 
application of statistics to the analysis of computational error. R. Creighton Buck, in “Linear Spaces 
and Approximation Theory,”’ discusses approximation theory from the point of view of functional 
analysis. In ‘‘Operational Methods in Numerical Analysis Based on Rational Approximation”’ Z. Kopal 
approximates non-polynomial difference and differential operators by ratios of polynomial operators 
and applies this to numerical quadrature and the numerical solution of ordinary differential equations. 
Philip J. Davis, ‘On the Numerical Integration of Periodic Analytic Functions,’’ compares the re- 
markably efficient trapezoidal rule with Gaussian quadrature and shows that widespread analyticity 
of the integrand favours the former, singularities near the integration interval the latter. Herbert E. 
Salzes states, in “Some New Divided Difference Algorithms for Two Variables,’’ criteria which such 
algorithms should fulfill and proposes two examples. Preston C. Hammer’s “Numerical Evaluation 
of Multiple Integrals’ is an expository paper on his and his collaborators’ recent results, especially 
on integration over spheres, cubes, cones and simplexes. The general problem studied in “Optimal 
Approximation and Error Bounds’’ by M. Golomb and H. F. Weinberger is stated by the authors to 
be “how to derive from a given amount of information about the unknown function u the total range of 
possible values of the functional F(u).’’ The authors discuss this general problem of approximation, 
particularly when and to what extent approximation is possible; they completely solve the case where 
the information consists of the values of a finite set of linear functionals /;(u) and a bound for a non- 
negative quadratic functional in uw. These results are practically evaluated and many special cases, 
applications and generalisations given. “The Rationale of Approximation’”’ by Arthur Sard is a general 
exposition of the theory of approximation of functionals. The following problem in Tchebycheff approxi- 
mation is, for instance, considered: how to choose a linear operator 7’ from a given linear manifold of 
bounded linear operators on a Hilbert space X such that ||7'g — h|| attains a minimum over 7’ for fixed 
g and he X. Joseph L. Walsh writes “On Extremal Approximation”’ and reviews recent results on the 
structure of extremal polynomials. In “Numerical Methods of Tchebycheff Approximation” E. Stiefel 
first outlines the basic theoretical tools required and then discusses old and new algorithms for solving 
the classical problem of linear Tchebycheff approximation. A survey of recent interpolation methods 
available to the user of tables and desk machines is contained in L. Fox’s ““Minimax methods in table 
construction.’’ T. S. Motzkin discusses very briefly and without proofs the ‘‘Existence of Essentially 
Nonlinear Families Suitable for Oscillatory Approximation,”’ in particular of families which have some 
of the advantages of polynomials for approximating given functions, ‘‘On Variation Diminishing Approxi- 
mation Methods,’ by I. J. Schoenberg, recapitulates some of the author’s work over past years. Let 
g(x) = L(f(x)) be a method furnishing the approximation g(z) to f(x). If v(f) is, in a certain sense, the 
number of changes of sign of f(z) in the range [0, 1], then the method is said to be variation diminishing 
if it has the property v(g) < v(f). The author shows, for example, that approximation by Bernstein 
polynomials has this property and thus avoids unnecessary oscillations of the approximating function 
such as often appears in least squares and Tchebycheff approximation; he also investigates approximation 
of a given function by power and by trigonometric polynomials from this point of view when the approxi- 
mation is generated by a linear functional operator. M. Golomb, in ‘‘Approximation by Functions of 
Fewer Variables,’’ investigates both least squares and Tchebycheff approximation of a function of n 
independent variables by functions of fewer variables, thus attacking a difficult practical problem and 
giving hope for the emergence of efficient numerical procedures. Practical applications of polynomial 
approximation are dealt with by J. C. P. Miller in “Extremal approximation—A Summary.” 


(Continued on p. 334) 
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A TWO-POINT METHOD FOR THE NUMERICAL SOLUTION 
OF SYSTEMS OF SIMULTANEOUS EQUATIONS* 


BY 
W. M. KINCAID 
University of Michigan 


1. Introduction. This paper is concerned with the problem, which arises frequently 
in practice, of finding one or more solutions of a set of n rig mana equations 
r,) = 0,7 = 1, --- , m, in m unknowns 2, , --- , x, . For convenience we 


Oitti s°** 93 
shall assume, as is usually true in actual cases, that the oui yg; are analytic in a 


region surrounding each solution of interest, although the approach we shall use would 
remain applicable under much milder conditions. 

While a wide choice of methods is available when the equations are linear, the same 
cannot be said in the non-linear case. Simplifying the system by elimination is likely 
to be laborious, and may be impossible if the equations are not algebraic. If the system 
is to be tackled in its original form, about the only general procedures are Newton’s 
method and the method of steepest descent. Both of these require the repeated evalua- 
tion of all n° partial derivatives 0¢,;/dx,; . In addition, Newton’s method involves the 
solution of a set of linear equations at oil step. The method of steepest descent does 
not require this, but generally exhibits only first-order convergence, whereas the con- 
vergence of Newton’s method is of the second order in general. 

These remarks do not mean that existing methods are useless, but they make it 
clear that there is room for fresh suggestions. In the following section we shall presen- 
an iterative method, believed to be new, that does not require the calculation of derivt 
atives or the solution of sets of linear equations, and yet displays second-order con- 
vergence. The method may be regarded as a generalization of the classical method 
of false position for the solution of one equation in one unknown. 

While experience with the method has not been extensive, it has been sufficiently 
encouraging to suggest that the method is worthy of serious examination. It may even 
be worth considering for the solution of linear systems. In this case the solution is 
attained (apart from rounding-off errors) after a finite number of steps, and the volume 
of computation appears comparable in magnitude with that required by existing methods. 

2. Description of the method. Since the method of solution to be described is a 
generalization of the method of false position, as stated earlier, a brief review of the 
latter is in order. 

Let f(x) = 0 be an equation in one unknown, and consider first the case where f 
is linear. If the value of f(x) is computed for two distinct values a and b of x, its value 
for all x is given by the identity 


f(z) = ~*~ [( — af() + (b — a)f(a)). 


*Received September 9, 1959. This work was conducted under contract to the United States Army 
Signal Air Defense Engineering Agency, Contract No. DA-36-039 SC-64627 
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Setting the right member equal to zero, one obtains 


_ af(b) — bf(a) 
f(b) — {@ 


for the solution of the equation. 

If f is non-linear, one can perform the same operation, obtaining the zero of the 
linear function coinciding with f for x = a and x = b. The operation can then be repeated, 
replacing a or } by the value of x thus obtained. Continuing in the same manner, one 
obtains a sequence of values of z that will converge to the desired solution if the function 
is well behaved and the starting points a and b are suitably chosen. The order of con- 
vergence of this procedure, which is probably the most efficient form of the classical 
method of false position, has been shown to be (1/2) [1 + (5)'””] (see [1, 2]) and is thus 
greater than 1. (By the order of convergence is meant the greatest value of k for which 
the ratio 6,/6:_, is bounded, where 6,-, and 6, are the departures of two successive 
iterates from the true solution). 

Now let f be a real-valued function defined on a region of Euclidean n-space in- 
cluding the points A(a, , --- , a,) and B(b, , --- , b,). By analogy with the one-dimen- 
sional case, we may determine a new point X having the coordinates 


a,f(B) — b,f(A) 
f(B) — f(A) 


(provided f(A) # f(B), as we shall assume in similar cases henceforth). For convenience 
(1) may be written in the abbreviated form 


X = AfB. (1a) 


It is apparent that the point X thus determined lies on the line joining A and B, 
and it follows from a slight generalization of our earlier statements that f(X) = 0 if f 
is linear; speaking geometrically, X is then the intersection of the line AB and the hyper- 
plane f = 0. Moreover, if A and B both satisfy another linear equation g = 0, the same 
is true of X. 

Making use of these ideas, one can set up a variety of methods for solving sets of 
simultaneous linear equations. One such method has been proposed by Purcell [3]. 
Some of the methods may be applied to sets of non-linear equations as well. It is appro- 
priate to refer to these as ¢wo-point methods in view of the character of the operation 
defined by (1) on which the methods are based. The remainder of the present paper is 
concerned with one method of this type, which possesses some particularly desirable 
properties. 

In the following we shall confine ourselves to the case of two equations 





xX; = (¢ = 1,2, --- ,n) (1) 


,(X) a o (2, » Xa) = 0, t = i. 2, 


in two unknowns, as it is illustrative of the general case. Without loss of generality, 
we may suppose that O(0, 0) is a solution. We shall assume that the vectors grad ¢,(0) 
and grad ¢.(O) are non-zero and non-collinear, so that the solution is a simple one. 

We begin the process of solution by selecting three linear combinations f{(X), g(X), 
h(X) of ¢,(X) and ¢.(X) in such a way that {/(X) + g(X) + A(X) = 0, while any pair 
of f, g, h are linearly independent. Thus the equations f(X) = 0, g(X) = 0, h(X) = 0 
have the common solution 0, and the vectors grad f(O), grad g(O) and grad h(O) are 
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non-zero and point in three different directions. The reason for introducing an additional 


equation will be indicated in a moment. 
We select three initial points R, , S, , 7, . Making use of the operation defined by (1), 


we locate successively the points 
Si = R,fS, , T; = R,fT, , R; = Sigh, , (2) 
Z3 = S{gT! ; Re = Th : ; S2 = ThS{ ° 


The equations (2) represent one iterative cycle. The next cycle proceeds in the same 
way, with R, , S, , T, in place of R, , 8,,7,. 

If the functions f, g, A are linear, the lines S{7{ and R,T, coincide with f = 0 and 
q = 0O respectively; thus 7, (also R, , S,) coincides with O. If the functions are non- 
linear, but not too badly so, one expects these relations to hold approximately, and 
in this case R,S, should coincide approximately with h = 0. Thus if all goes well, suc- 
cessive cycles should yield a sequence of approximately similar triangles converging on O. 

The fact that the sides of the triangles approach distinct fixed directions is important. 
For otherwise one of the triangles might collapse to a line, in which case all succeeding 
points would lie on the same line, and the solution could not in general be attained. 
It was with these considerations in mind that the two original equations were replaced 


by three. 
The foregoing discussion is of course quite heuristic. We shall show later that it 


can be made rigorous if the initial points satisfy suitable conditions. 
3. Useful definitions and lemmas. Interpreting a point P in the plane as a vector 
from O to P, we may define the linear functions 


u(P) = P-U, vo(P) = P-V, w(P) = P-W, (3) 
where 
U = grad f(0), V = grad g(0), W = grad h(O). (4) 
Thus u(P), v(P), w(P) coincide respectively with the linear terms of the Taylor ex- 
pansions of f, g, h about O. Hence 
u(P) + o(P) + w(P) = 0. (5) 
We see that any pair of the three numbers u(P), v(P), w(P) may be regarded as a linearly 
transformed set of coordinates of the point P; this interpretation will underlie our 
subsequent discussion. The ‘axes’ u = 0, v = 0, w = O are the respective tangents to 
the curves f = 0,g = 0,h = Oat O. 
The non-linear parts of f, g, h are given by the functions 
F(P) = f(P) -— uP), GP)= 9P)-P), HP) =P) — wl). (6) 
For any two points P, Q in an appropriate neighborhood of O, we have by the mean 
value theorem 
F(P) — FQ) = (P — Q)-grad F(6P + (1 — 4)Q) (7) 
for some @ in [0, 1]. Since grad F(O) = 0 and F is analytic (in z, y and thus in u, v) 
near O, there exist positive constants A’ and K’ such that 


| grad F(P) | < A’ max (| u(P) |, | o(P) |) (8) 
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whenever max (| u(P) |, | v(P) |) < K’. It follows from (7) and (8) that 
F(P) — FQ) | < A’o(P, Q) max (| u(P) |, | u(Q) |, | oP) |, | x) J), (9) 


(where p(P, Q) is the ordinary Euclidean distance between the points P and Q) for 
max (| u(P) |, | w(Q) |, | o(P) |, | 0(Q) |) < K’. 


These relations make it convenient to define the norm || P — Q || of a vector P — Q 
by the relation 
'P —Q|| = max (| u(P) — u(Q) |, | oP) — o(Q) |, | w(P) — w(Q) })- (10) 
In particular, we write 
P || = ||P -—O]| = max (| u(P) |, | oP) |, | w) )). (11) 
Clearly the norm, so defined, has the usual required properties, and the ratio 


|| P — Q||/p(P, Q) has positive upper and lower bounds. Thus from (9), (10), and (11) 
we conclude that there exist positive constants A and K such that 


F(P) — F(Q) | < A || P — Q || max (|| P ||, || Q |) (12) 
for max (|| P ||, || Q ||) < K. In particular, 
FREI S AUP fo [Pil < &. (13) 


tepeating the same arguments, we may suppose A and K so chosen that (12) and (13) 
hold with G or H in place of F. 
We note at this point that (5) and (11) can be combined to yield 


P || = min (| w(P) | + | oP) |, | uP) | + | w(P) J, | oP) | + | we) J); (14) 


a similar relation holds for || P — Q 
Finally, we introduce 
u(P) — u(Q)| | uP) — u(Q) ) is 
| r ”) 


ri(P,Q) = max ( —A— 
iis v(P) — v(Q) |’ | w(P) — w(Q) | 


with similar definitions for r,(P, Q) and r,,(P, Q). We see that r,(P, Q) may be regarded 
as measure of the difference in direction between the line PQ and the line u = 0; paral- 
lelism corresponds to r,(P, Q) = 0. 

From (10) and (15), we have 


|} uP) — u(Q) | <r, Q) ||P — QI. (16) 
On the other hand, if »(P) # v(Q), one obtains from (14) and (15) 


IP —Q|| < | uP) — u(Q) | + | oP) — ro(Q) 


| u(P) — u(Q) 
= | oP) — °(Q) | (| — — ee 7 
a 7 ( v(P) — v(Q) sl (17) 


l 
| 


< | o(P) — o(Q) | (P,Q) + 1). 
The inequalities (16) and (17) are clearly valid for any permutation of the letters u, 
v, w. Geometrically, (17) implies that if PQ is nearly parallel to u = 0, the ‘distance’ 
between P and Q can be only slightly greater than the difference between their v- (or w-) 


coordinates. 
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4. Convergence of the process. In view of the foregoing remarks, the sequence 
of points defined by (2) will converge to the desired solution if and only if 


lim || 2, || = lim {| S, || = lim || 7, || = 0. (18) 


no n—-+@ 


We shall prove the following result. 
Theorem. Under the assumptions stated previously, the sequence defined by (2) 
will converge to the solution at O if the initial points R, , S, , 7, , satisfy the conditions 


7 1 
ax (IIR. WS. I 1 II _ 3K ( 
max (|! R, ||, || Si {], 1] 7, |) < min ee 4 ) 4 (19) 


where A and K are the constants appearing in (12), and 


max [r.(S, , 7), 7.(Ri , Tr), Po(Ry , S)] < TT (20) 


Geometrically, the conditions (19) and (20) mean that the points R, , S, , 7, lie 
within a specified neighborhood of O and that the sides of the triangle R,S,7, are approxi- 
mately parallel to the lines u = 0, v = 0, w = 0. As will become evident, other com- 
binations of constants could be used in (19) and (20); the important point is that some 
sufficient set of bounds can be given. 

Writing 7, max (|| R, ||, || S, |], || 7. ||) forn = 1, 2, --- , we see that to estab- 
lish the theorem it is sufficient to show that the stated hypotheses imply (a) that there 
is a constant p < 1 such that M, < pM, , and (b) that (20) is satisfied with R, , S, , 7; 
replaced by R, , S. , 7, . For if (a) and (b) hold, the same argument can be used to 
show that 1/7, < pl/,,M, < pM, , ete., from which (18) follows immediately. 

As one would expect, the proof requires repeated use of the inequalities (12), (16), 
and (17). For (12) to be applicable, the points involved must have norms <K, as is 
true of the initial points R, , S, , 7’, by hypothesis. For (16) and (17) to be applicable, 
directional conditions analogous to (20) must be satisfied. Defining 


M,, = max (|| R, HU SPL, UT Ip, (21) 
M,. = max (|| Ri ||, || Sf |!, || 7. I), 
we see that these requirements may be met by establishing the inequalities 


0 | Ma<K, rv (RT, <i. (22) 


M,, <K, r({S; ,Ti) < 10 


An outline of the proof will now be presented, followed by a more detailed discussion 


of some specific points. 
Making use of the properties of the initial points, we begin by deriving the results 
21AM, 


ital nl!) ' | w\ | 1] > | Siena enti “ i | | 
max [| w(S7) |, | u(7) |] < || RB, | i0 — 114M, < 0.111 || R, |], (23) 


and 


| 5 ies ll 214M 
ax [| w(S%) |, | (7) |] < j=+ 1 | 
max [| w(Si) |, | o(7%) |} S I Rs II E 10710 — Moan | SURI, A) 
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from which follow, in view of (14), the relations 


M,, < 1.22 || R, || < 1.22M, < .92K, AM, < .0611. (25) 
The next step is to establish the inequality 
11AM 
. S! N< PR kek . q- ; %6 
r(S{iT)) < 9— 114M, < .0651 (26) 


The results (25) and (26) permit the relations (12), (16), and (17) to be applied 
to the triangle R,S{7{ , in line with the reasoning outlined earlier. 
In an analogous manner, one obtains the further inequalities 


My. < 1.251 || R, || < .94K, AM,, < .0626, (27) 

and 
, 11AM, “oe ” 
r(RiT.) < 9 — 1AM, < .0807, (28) 


thus showing that the triangle R{,S{7; also has the desired properties. 
Since 7, would coincide with the solution point if the functions f and g were linear, 
one would expect || 7. || to be sharply bounded. In fact, one can show that 


r(Si71) || Ti |] 1 + A | TF II) 


1 — AM,,[1 + r.(S;7))] 


: .0651(1 + AM,,) . 
| iis | ‘ 
111 {Ra ll + 1 TE | opera S 201 1 Ri Il (29) 














| u(T.) | < | (7) | + 





lA 


and 


| o(Ts) | < || 82 |] tA << 1168 || R, || (30 
| VV 277i —_ i 1 10 — 11AM,, —° it 1 {i . ) 


so that by (14) we have 
[| T2 || < .369 || R, ||. (31) 
Identical arguments, applied to the triangle R,S.T, , yield the inequalities 
| Ro || S .254 || Ri |l, |] S2 || < .282 || R, |I, (32) 
which together with (31) imply the relation 
M, < .369M,. (33) 


Thus assertion (a), stated at the outset, is verified. 
By a further repetition of previous arguments, one obtains 


11AM,, 
a | ( d 
9 114M, < 0829- (34) 





To(R2S2) S 


Clearly (26), (28), and (34) constitute a verification of assertion (b). Thus the proof 
is completed. 

To a large extent, the arguments used in establishing successive steps of the preceding 
outline are variations on a single theme. Accordingly, only one or two steps will be 
considered in detail. 
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To begin with the proof of (23), we have in view of (2), (3), and (6) 








(st) = YedKS) = u(S,)f(Ri) __ __u(R) F(S:) — u(S,) F(R) 
~ f(S,) — f(R,) u(S,) — u(R,) + FCS, — — 
7 COE , 





1+ 2S) — wR) 
3y (12), (17), (19), (20), and (21) we obtain 














|F(S:) ) — F(R,)| _ F(S,) — F(R) | _ || 8, — RII 

u(S;) — u(R,)| |S, —R, |] | uGS,) — u(R,) | 
< ie (|| R, II, I] S, ||) }f2 + r.(R,S,)] (35) 
<i m AM, . 

Combining these results with (11), (13), and (19) gives us 

ju) | AM, + A || RIP en 
| u(S;) | < ——_! < Il Ri Il p= 1AM < 0.111 || R, |], 
1 — ie AM, ie 


and repeating the argument with S{ replaced by TY yields (23). 
The proof of (24) starts with the decomposition 


w(R,) — w(S,) 
u(S,) — u(R;) 





(u(R,) + F(R,)] 





w(S{) = w(R,) + F(S,) — F(R.) , (36) 
u(R,) — u(S;) 


and proceeds in much the same way. The proof of (26) likewise involves similar argu- 
ments, starting with the identity 
r | 
|u(R,) — u(T,) f(R,) — f(T.) | 


u(S) — u(T}) _ | u(S,) - u(T,) f(S) — (1) | - 
S}) — o(T)) o(R,) — o(7') F(R) — f(T) 
u(S,) — o(7,) f(S:) — f(T) 

The proofs of (27) and (28) require the use of (25) and (26), but involve no new 
ideas. The proof of (29) begins with a step analogous to (36), and uses the bound on 
| u(T) | given by (23). The remaining steps fall into the same pattern. 

Corollary. Under the hypotheses of the preceding theorem, convergence is of the 


second order at least. 
This result is easily established. Referring to (23) and replacing AM, by its bound 


0.05 in the denominator, we obtain 


| u(S)) | < 2.224M? ; 








likewise, 
| u(Ts) | < 2.22AM;. 
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Similarly, from (26) 
r(SiT)) < 1.30AM, . 
Substituting these results in (29) and taking account of (25) and (26) yields 


1.30AM, (1.0611) 


m\ } ‘ } 2 9 PP tn tote Rowe, SL \ that tet ak até 
w(T's) | S$ 2.22AM; + 1.22M,-7 7 o¢57 (611) 


4.02AM; . 


Il 


Since from (39) 
| o(7.) | < 2.17AM; , 
we get 


| 7. || < 6.19AM; . 
Similarly, 
| R. || < 3.61AM; , || Se || < 4.48AM?. 


Thus M, < 6.19AM; , which establishes the corollary. 

Actually the convergence appears to be slightly faster than second-order, since 
w(R,) and w(S,) are of the order of M? . 

5. Numerical example. To illustrate the application of the foregoing ideas, we 
may consider the set of equations 


f(x, y) = x’ — 4y = 0, gz,y) =y — 2x+ 4y = 0, (38) 


to which we may add 
h(x, y) = —f(z, y) — g(x,y) = —x* —y +22 = 0. (39) 


One of the common solutions of these equations is (0, 0). Taking for starting points 
R,(0, 1), 8:01, —2), 7:(—1, —1), and applying (2) yields the sequence of points listed 
in the table and shown on Figs. 1, 2, and 3. The accelerating character of the convergence 
is evident, as is the tendency of the points to form similar triangles. 

One aspect of the computation deserves further discussion. The point 2, , which 
is conceptually an approximation to the point of intersection of the line R/T, and the 
curve h = 0, happens to fall very close to the curve f = 0. In consequence the points 
R, , S} , T are much closer to each other than to the solution point, and a considerable 
extrapolation occurs in determining the points Rj and 7; . While no difficulty arises 
in this example, it is apparent that if R, had actually fallen on f = 0, the points S} 
and 7’; would have coincided with R, , and R{ and 7, would have been undetermined; 
in a practical computation, where round-off errors must be allowed for, something 
less that exact coincidence could cause serious trouble. A number of expedients may 
be imagined for meeting this problem, and are now under consideration. 

From (38) and (39) we find readily 


u= —4y, vy = —2xr + 4y, w = 22, (40) 


and 
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TABLE 
Example of Iterative Solution 
1 y ys q h 
T, —1.0 —1.0 5.0 
Si 1.0 —2.0 9.0 
R, 0.0 1.0 —4.0 5.0 
Ti — .444 lll 1.344 
St .308 .077 — .302 515 
Ri 290 .130 .479 
Ts 170 .083 — .303 .304 
S2 0290 0917 — .366 
R, — .0389 .0012 — .00329 .0826 
T’ — 0412 .00030 .0836 
FA - 0390 .00038 0795 — .0795 
Ri — .0413 — .0206 — .0847 
‘ig 00388 00194 — .00774 .00774 
S; 0,760 .00180 — .00720 
R; 0,978 .0,497 — .0;199 
aT - .0;190 — .0;,97 
S3 0,984 — .0;23 








Since for any values of w, and w, we have 
|4w, —3w,| = 3 | wv, —w.||w, +w.! < $+] wv, — w, | max (| w, |, | w, |), 
we conclude, referring to (10) and (11), that for any two points P, Q 
F(P) — FQ) | < 3 ||P — Q || max (|| P [], |! Q/)). 
Similarly, we find 


G(P) — G(Q) | 


lA 
lm 


|P — Q|| max (|| P |i, || Q])), 
and 
H(P) — H(Q)| < § ||} P — Q]! max (|| P ||, || Q])). 
Comparing these statements with (12), we see that the requirement (19) becomes 


max (j] R, ||, || Si |], 1) TID) <5, 


) 


or, if x, y are the coordinates of any one of the three points R, , 8S, , 7, , 


! 1 1 2 
max (| 27 |, | 4y |, | 2a — 4y |) < On. 

Substituting from (40) into (15), we can place limits on the directions of the lines 
R,S, , R,T, , and 8,7; , corresponding to (20); for example, we find that the bounds 
on the slope of 8,7’, are —1/20 and 1/21. 

It is apparent that the bounds suggested by the theorem proved earlier are much 
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too conservative in this case. However, the result does show that the process will con- 
tinue to converge once the solution has been closely approximated, and indicates the 
character of the sequence of iterates. 

6. Concluding remarks. The foregoing example indicates that convergence may 
take place in practice under much less stringent conditions than were postulated in 
the statement of the theorem. Apart from the possibility of a more refined analysis, 
one must recall that in the course of the proof we necessarily assumed that each quantity 
appearing took on its most unfavorable value at each step, which would be highly 
improbable in an actual example. On the other hand, it appears that some conditions 
of the type imposed are necessary for convergence to the desired solution. For if the 
location of the initial points were uncontrolled, the iterative process might converge 
to another solution of the system, while their orientation must be restricted to keep 
the denominators in (2) away from zero. 
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The iterative process :aay clearly be generalized to apply to systems of n equations; 
presumably the same is irue of the convergence properties. Moreover, we may note 
that the foregoing discussion is applicable to the finding of complex as well as real 
solutions. 

As remarked in Sec. 2, the method of solution discussed in this paper is only one 
of a variety of two-point procedures that may be constructed. Some of these have been 
tested successfully on examples, but no sufficient conditions for convergence have been 


established. 
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CONDUCTION-CONVECTION FROM A CYLINDRICAL SOURCE 
WITH INCREASING RADIUS* 


BY 
H. R. BAILEY 
The Ohio Oil Company, Littleton, Colorado 


Summary. The problem of heat flow by conduction and convection from a 
cylindrical source with increasing radius is solved. A quasi stationary state solution is 
obtained for the case of a finite convection coefficient and with the radius increasing at 
a constant velocity. A transient solution is obtained for the case of an infinite convection 
coefficient and with the radius increasing at a rate proportional to the square root of 
time. 

In the latter case an explicit evaluation of an integral form of the solution is obtained 
by showing that the solution, in a certain coordinate system, is independent of time 
and thus the partial differential equation reduces to an ordinary differential equation 
which is solved explicitly. 

Introduction. The theory of heat flow from a moving source is of interest in a 
number of applications; some of these are discussed by Crank [5]t. Most of these appli- 
cations are concerned with heat conduction from a moving plane source. 

The problem of heat flow from a cylindrical source with increasing radius has been 
studied by H. R. Bailey, B. K. Larkin and H. Ramey, [1], [2] and [9]. This problem arises 
in connection with a secondary oil recovery process by underground combustion. The 
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Figure 1 





*Received September 16, 1959. 
t}Numbers in brackets refer to the references listed at the end of the paper. 
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papers mentioned above are mostly concerned with the conduction mechanism for heat 
flow; however, Ramey [9] does give an approximation of the convection effects. Heat 
flow by conduction and convection, without moving sources, has been studied by Bland 


[4] in connection with heat exchanger applications. 
In the present paper we consider heat flow by conduction and convection from a 


cylindrical heat source with increasing radius, r = r; , as shown in Fig. 1; gas (air) is 
injected at r = 0 into a porous solid of infinite extent and it is the oxygen from this 


air that supports the combustion. With certain simplifying assumptions (see Bailey and 
Larkin [3]) the equations describing this process can be written in the form 


) aT . 
+k’, d =k'1T, —T) = PX) aT : (1) 


a, Lan, at, 
r or 


or” r or ot 
where 7’, is the temperature of the solid, 7’ is the temperature of the gas, k is the con- 
ductivity of the solid, a’ is the reciprocal diffusivity, h is the convection coefficient, r is 
the radius, ¢ is time, ®(r, ¢) is a source function and P(t) is a function of time whose 
form depends on the air injection program. From the second equality in system (1), we 
have 


7 


kP(t) aT 


_ 7 a wt 
T= T+ = T+ EG, OS 


- 
where the identity is a definition of E(r, t). With 7, replaced by the above expression, 
system (1) can be written as a single partial differential equation: 


a(T + EaT/ar) , 1 A(T + E aT/ar) 


_- = 


2 | 
or r or P 
(2) 


Z cry 7 4/1 /4,.) : 7 ; IT 
— g@ METS 6) + bag, 9 — ene = 0. 
ot or 








The functions P(t) and ®(r, ¢t) depend on the assumed air injection program; two 
cases are of particular interest in underground combustion: (a) constant velocity, 
r,; = v,l, where v, is a constant, in this case it is shown in [3] that P(t) = k, vit and 
&(r, 1) = qv,d(r — r,), (b) constant injection rate, r; = 2Ut, where U is a constant, in 
this case it is shown in [3] that P(t) = 2n where n is a constant and ®(r, t) = qUr;' 
5(r — r,). In the above equation k, and q are constants and 6(r — r,) is the Dirac delta 
function. 

In Sec. 2 of this paper a quasi-stationary state solution of Eq. (2) is obtained for the 
constant velocity case. Solutions of Eq. (1) are obtained in Sec. 3 for an instantaneous 
cylindrical source and for a cylindrical source with increasing radius (i.e., the constant 
injection rate case). This latter solution is obtained in the form of an integral which is 
evaluated explicitly in Sec. 4. 

2. Quasi-stationary state solution for the constant velocity case. In this section 
a quasi-stationary state solution of Eq. (2) is obtained for the constant velocity case, 
i.e., 7, = vst with v, a positive constant. Clearly there would not be a steady state solution 
in the usual sense; however, it is assumed that an observer moving with the source 
would see steady state attained for sufficiently large values of time. This assumption is 
discussed by Jacob [7]. 
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Putting r = s + v,t transforms Eq. (2) into 
2 7 Al ™/5 e Dy y 
o(T + E d77/ds) 4 1 0o(7 + E dT/ds) + oy, o(T + E dT /ds) 














as" s + ot ds 08 (3) 
20T +E AT/ds) hE OT , yy. _ 


where E = k P(t)/h(s + v,t) = kk, vjt/h(s + vf). The variable s measures the radial 
distance from the moving source, and the assumption that quasi-stationary state is 
attained means that 7'(s, t) does not change with respect to time for sufficiently large 
and finite s. Thus, passing to the limit as t + ©, Eq. (3) is reduced to the ordinary 
differential equation, 


ptt , eT 


ai 2m, T 
* ds® +r ds* + ake, 


2 
a +av 
ds” . 





a — ky, a = 0, (4) 
where the last term in Eq. (3) is replaced by the equivalent boundary condition, (a), 
below. The term k,v, dT’/ds is the limiting form of the corresponding term hE/k dT/ds 
in Eq. (3), this follows from the definitions in Sec. 1 of E(r, t) and P(#) for the constant 
velocity case. 

The boundary conditions are: 














) @:| _at.|. _aT| aT 
. d8 |,<04 Sine Glu G8 \eae- 
. aT = aT a —Qvy 
+& rina fn 6 * 
> | » | 
(b) dT | os dT 
2 OF tien 

(c) T remains bounded ass— — o~, 

(d) T—Oass— — o, 
where dT/ds |,.,, and dT/ds |,.o- are the right and left hand derivatives respectively 
ate = GO. 


Condition (a) is equivalent to the source term in Eq. (3). Condition (b) results from 
the assumption that d7/ds is continuous at s = 0, i.e., that the heat source is only in 
the solid phase. Conditions (c) and (d) follow immediately from heat balance considera- 
tions for a moving source. 

The general solution of Eq. (4) is of the form T = Cy + C, exp (r,s) + C2 exp (res) 
and it can be seen from the characteristic equation determining r, and r, that they are 
both real and negative, provided /, is less than a’ which is the case in underground 
combustion. Thus, in order to satisfy boundary condition (c) for s < 0, we must have 
C, = C, = 0. For s > 0, and to satisfy boundary condition (d), we must have Cy = 0. 
Imposing boundary conditions (a) and (b) results in the following solution of (4) subject 
to conditions (a), (b), (ce) and (d). 





Rq 
T=——t., s<0 
(R — S)a’k (5) 
2Rq {2 [((-1 — R + W)s/2E] _ exp [(-1 —R — Wie/2E I} eT 
~ @kW i+R—-W 1+R+W ne 


where R = a’Ev,, S = k,Ev, W = [(1 — R)? + 48]'”. 











328 H. R. BAILEY (Vol. XVIII, No. 4 


@ 


TEMPERATURE FRACTION T 
ES o 





0 | 2 3 4 5 6 7 8 9 10 
DIMENSIONLESS DISTANCE AHEAD OF FRONT s\;/a 


Figure 2. 


The above solution is also a solution for the corresponding linear problem of a plane 
source moving at a constant velocity in a semi-infinite medium. This follows from the 
observation that a cylindrical source approaches a plane source for sufficiently large 
radii. 

If it is assumed that conduction is the only mechanism for heat transfer then the 
corresponding formula for the temperature has been obtained [1] by determining the 
limit as t — © of an explicit solution of the partial differential equation. The latter 
method is more difficult; however, it is more satisfying since it is not necessary to make 
the a priori assumption that a quasi-stationary state will be obtained. 

It is shown in [3] that R = 5.18 S if the injection gas is air; Fig. 2 shows the temper- 
ature profiles for this case. The temperatures are given in terms of the temperature 
fraction, r = T'/(q/a’k). 

3. Transient solutions for the constant injection rate case. In this section explicit 
solutions of Eq. (2) are obtained for the constant injection rate case, r; = 2Ut, assuming 
infinite convection coefficient, h. With these additional assumptions we have EF = 0, 
and P(t) = 2n, where n > 0 is a constant depending on the parameters in the physical 
problem. Equation (2) thus takes the form 

OT ,1—20T .,aT 


ar a a —a Ot +k P(r, t) = (0. (6) 


In the following analysis two forms of the source function, ®, are considered: (1) an 
instantaneous cylindrical source, (2) a continuous cylindrical source with increasing 


radius. 
3.1. Instantaneous cylindrical source. Consider a cylindrical source with radius ro 


i 
| 
| 
i 
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which liberates one unit of heat per unit area at time ¢ . In this case ® can be represented 
in terms of Dirac delta functions and Eq. (6) becomes 
er, 1—2naP 
or” r or 


a 2 +k 8 — 1) Ht — &) = 0. ” 


The substitution 7 = Sr” transforms the above equation into 
a’s ] as n? 2 os 


soap ire ane == «= 
or” r or r* ot 


+ k'r'™ 6 — ro) H(t — to) = O. 
Let R be the Hankel transform of S, e.g., [10], i.e., 


R& = [ rJ,.(ér)S(r, t) dr. 


Then taking the Hankel transform of the above equation gives P 
?R — a’ a = —k"' [ r'-"J,(ér) &(r — ro) &(t — to) dr 
“0 


= —k rs" J (eto) 6(t —_ to). 
Solving the above equation we obtain the particular solution 
R = k'a~*r, "J (to) exp [#(t — t)/a"] 


and taking the inverse Hankel transform, e.g., [6], gives 
Sr, t) = | ERE, DJG) dé 
= 2°k'(t — to) ro” exp [—a°(r? + 13)/4(t — to), [rroa”/2(t — t)]. 


Finally, from the transformation relating S and T, the following formula for the temper- 
ature is obtained: 
Tie = 20k '(t — te) "ror" exp [—a°(r? + 15)/4(t — to) J,[rroa’/2(t — t)], (8) 


where the subscripts ic have been added to T to indicate temperature due to an instan- 
taneous cylindrical source. 
For n = 0 the above expression for 7’;, reduces to the known, [1], solution for con- 
duction only. For n > 0, it is indicated below that the solution has the required properties. 
It can be shown by direct computation that 7;, , for t > 0, satisfies 
aT 1 — 2n oT _ 2 or 


or” r or i Sa 





0. 


Also, from the asymptotic form of J,(z) for large z, it can be seen that 7;. — 0 when 
t — t, at all points except on the cylinder r = ro where 7';, becomes infinite. By using a 
generalization of Webber’s first integral [11], 
ace b” 
2 2 n+1 
xp (— all du = = 
I exp (—a'u’)J,(bu)u (2a) 


XN 





i exp (—b*/4a’), 


we have for it > Zt, 


a’k [ QarT ;.(r, t; To , to) dr = Qo 
0 


which is the amount of heat liberated per unit length by the instantaneous source. 
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3.2. Cylindrical source with increasing radius. It is reasonable, [2, 3], to assume in an 
underground combustion process that the source function can be expressed by the 
formula 

@(r, t) = g or — r,) dr, /dt, 
where r,(¢) is the radius of the cylindrical source, 6 is the Dirac delta function and q 
is constant. The case of most practical interest is with constant air injection rate; and 
in this case the position of the front is given by the formula, r7 = 2Ut, where U is a 


constant. 
With the source function described above Eq. (6) becomes 





aT _1—2o0T 207 , qU 

—— —_————— —— — 9° — +4 i(ir—r,) = q 9) 
or” T r or "a ° & a v1) 0 , 

The solution of the above equation with 7’ = 0 at ¢ = 0 is given in terms of 7;, , Eq. 


(8), by the formula 


at 
T(r, t) = qU | [ T(r, t3 70 , tor; Or — r,) dro dbo , 


70 


where 7, is evaluated at t) , i.e., r, = (2Ut,)'’*. Performing the indicated integration 


with respect to rp we obtain 

T = 2°k'qUr" | (2Ut)-""(t — &)~ 
70 (10) 
-exp [—a°(r? + 2Ut))/4(t — t)]I,[27'\(t — to) *a7r(2Ut,)'”?] dtp. 


4. An explicit evaluation of the integral in Eq. (10). An explicit evaluation of the 


above integral, Eq. (10), is obtained by applying a method introduced in [1] for the 
case of conduction only. Making the substitutions y = r/r,; and r = (t — t)/t, Eq. 


(10) is transformed into 


1 
T(y, 1) = 2°"k"qUy" [ 7 (1 — 7)” exp [—@’'Uly’ + 1 — 7)/27) (11) 
70 
[7 a?Uy(l — 7)'”’] dr. 


Thus 7'(y, ¢) is independent of ¢ since ¢ does not appear in the right side of the above 


equation. 
The transformation y = r/r; in the partial differential equation (6) gives 
oT pant ao 2 oT 
os ty — — %°U - 
~ roe, ay 1 2 UNS 2a = 0, 


where the source term is replaced by boundary condition (b) below. Since we have shown 
that 7(y, ¢) is independent of ¢t we have d7'(y, t)/dt = 0 and the above equation becomes 


aT 1 — 2n ar dl 
9 
dy? : “he 2By ‘> = 0 (12) 





where B = a’U/2. 
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The boundary conditions are 
(a) =T, at y=0 


ar| _ a 
dy |yeis dy 








(b) = —qU/k 


y=l— 
(c:) T-0 as yoo. 

Condition (a) corresponds to the inlet temperature of the injected gas. It should be 
noted for the solution given by Eq. (11) that 7, = 0, and thus we have shown that 
T'(y, t) is independent of ¢ only for 7, = 0. We shall assume that 7'(y, ¢) is independent 
of ¢t for any value of 7’, and check the resulting answer to see that it satisfies the partial 
differential equation and the boundary conditions. Condition (b) replaces the source 
term, k~* @(r, t) = qUk™'r;* 6(r — r,); that is the heat source at r = r, is equivalent 
to a discontinuity in the temperature derivative at r = r,;. 

The ordinary differential Equation (12) can be solved either by making the sub- 
stitutions z = By’ and T = We-* which reduce Eq. (12) to a confluent form of the 
hypergeometric equation or by the substitution w = 07°/dz and integrating the resulting 
first order equation. The general solution of (12) is given by 


T C; + C;T(n, By’), y = 1 (13) 
T K2 + K;T(n, By’), y< 1 


where I'(a, b) is the incomplete gamma function defined by the integral 


I 


T(a, b) = [ e“u®* du. 
b 


Values of the incomplete gamma function have been tabulated, for example see Pearson 
[8]. The constants C, , C; , K, and K; are determined by the boundary conditions (a), 
(b) and (c), and the requirement that 7 be continuous at y = 1. These conditions lead 
to the following relations respectively: 
(a) T, = K, + K; T(n) where I'(n) is the gamma function 
(b) — 2C,e-7B" + 2K,e-"B" = — qU/K 
(continuity at y = 1) C,; + C; T'(n, B) = K, + K; T(n, B). 
Solving the above equations for C, , C; , K, and K; and combining with Eq. (13) 
we obtain 
e"B-"T(n, By” a es ae ae . 
Dn, BY) rp e-®B" 4+ 2k gU(T(n) — Tn, B)], y>1 
T'(n) (14) 


T = e®B-"qQU2"k"T(n, B) + Tee [T) — gU2"k'e*B"T(n, B)], y <1. 





T = 


It can be shown by direct computation that the above expressions for 7’ (with y replaced 
by r/r;, 7; = 2Ut) satisfy the partial differential equation (9) and the conditions T = 0 
att = 0, T = T, atr = Oand T—Oasr— o. 

Figure 3 has been obtained by evaluating Eq. (14) for the case T, = 0 and r/r, = 1. 
The fractional temperature rise, r = 7'/(q/a’K), is plotted as a function of a dimension- 
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Figure 3. 
less injection rate, 2B = a°U, for various values of a relative gas velocity y = n/B. 


The curve y = 0 corresponds to assuming conduction only. The value y = 0.2 is typical 
for an underground combustion process. 

Remarks. (1) For n = 0 the Equation (14) reduces to corresponding result, [1], for 
the conduction only case. In this case the terms in 7, do not appear, and thus the temper- 
ature at r = 0 cannot be prescribed but is determined by the solution. For n > 0 the 
temperature at r = 0 is prescribed and corresponds to the inlet temperature of the 
injected gas. 

(2) Since the solution is unique, we may equate the two solutions given by Eqs. 





(11) and (14) (for 7, = 0) and obtain the following evaluation 
[ em. <73 exp [—B(y? + 1)/7], (2uBO — 2") 
o (1 — 7) \ T 
: wit die. By’) —- Hie Bs | : y>1 
le-yre, B)| 1 ns ou BY ! y <i. 
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BOOK REVIEWS 


(Continued from p. $12) 


R. Creighton Buck presents a valuable ‘Survey of Recent Russian Literature on Approximation,’ 
with a bibliography of 135 items, most of them with their Mathematical Reviews reference. ‘The 
and Epsilon Algorithms’’ are compared by F. L. Bauer and formulas given con- 


Quotient-Difference 
continued fraction with its power series expansion. J. Barkley Rosser, in ‘Some 


necting the Stieljes 
Sufficient Conditions for the Existence of an Asymptotic Formula or an Asymptotic Expansion,” 
derives the asymptotic behaviour, for large z, of [2 f(x) exp (zg(x)) dz, giving numerical illustrations 
of his results. John W. Tukey, in ‘‘The Estimation of (Power) Spectra and Related Quantities’”’ presents 
various topics in the spectral analysis of time-series. Aliasing, resolvability and variability are discussed 
and the generalised Fourier analysis approach is contrasted with the ensemble approach. The salutary 
effect of prewhitening is described, as are various other computational techniques useful in the empirical 
determination of the spectrum. ‘‘Approximation in Partial Differential Equations,’”’ by L. Collatz, 
deals with the following problem. Let £ be a linear or quasi-linear second order differential operator 
of elliptic, parabolic or mixed (e.g. Tricomi) type, and D an n-dimensional region with boundary C. 
The problem is to solve E[u] = 0 in D with F[u] = 0 on C. It is attempted, in particular, to find a 
v satisfying E[v] = 0 only, such that maxp |u —v| is small, and examples of the technique are given. 
Finally, John Todd outlines a graduate course on the role of “Special Polynomials in Numerical Analy- 
sis,’’ embracing Weierstrass’ theorems, Tchebycheff polynomials, the theorems of the Markoffs and 
related topics. 

This volume will be greatly appreciated both by persons interested in the pure mathematical as- 
pects of approximation theory and by those interested in applying the theory to particular situations. 


WALTER FREIBERGER 


Théorie des Circuits Non-linéaires en Régime Alternatif. By V. Belevitch. Librairie Uni- 
versitaire, Louvain, Belgium, and Gauthier-Villars, Paris, 1959. 293 pp. $8.45. 


This is a sequel to an earlier text by the author, ‘Theory of Telecommunication Circuits,’’ which 
dealt primarily with linear time-invariant circuits of the type used in communication systems. The 
present text is concerned for the most part with the analysis of rectifiers, modulators and oscillators. 
The material is oriented toward the needs of the circuit designer rather than the mathematically in- 
clined engineer. There are hardly any differential equations in the text, the author’s main tools of analysis 
centering on the use of Fourier series. Although these are employed skillfully and effectively, the range 
of applicability of the author’s methods is quite limited. 

An interesting part of the book is a nine-page section describing in brief recent contributions to the 
subjects discussed in the text. One wishes that some of these contributions were discussed more fully 
in the main body of the work. 

All in all, this text is an impressive testimony to the author’s high competence. Its narrow scope, 
however, robs the book of much of its value to those who have a broad interest in the theory of nonlinear 


systems, 
L. A. ZADEH 


Methods based on the Wiener-Hopf technique. By B. Noble. Pergamon Press, New York, 

London, Paris, Los Angeles, 1958. x + 246 pp. $10.00. 

The Wiener-Hopf method is a fascinating and useful technique whose consistent exploitation in 
the treatment of boundary value problems and integral equations has evolved during the last fifteen 
years, despite its invention in 1931. The method depends only on classical theorems in the theory of 
functions of a complex variable and the author reviews all such necessary background material in 


Chapter I of this book. In Chapter II he utilizes three different formalisms, each depending on the 


(Continued on p. 354) 
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APPLICATION OF CONFORMAL MAPPING TO VISCOUS FLOW 
BETWEEN MOVING CIRCULAR CYLINDERS* 


BY 
LEE A. SEGEL 
National Physical Laboratory, Teddington, England 


Abstract. This work shows that conformal mapping provides an effective way to 
solve certain unsteady two-dimensional perturbation problems involving the flow of 
a viscous incompressible fluid, in particular flow between moving circular cylinders. 
If the outer cylinder is considered fixed, the principal motions treated are the slow 
rotation of a slightly eccentric inner cylinder, and the vibration of an inner cylinder 
about a slightly eccentric point. Mapping the given circular boundaries (of a cross- 
section) into concentric circles enables one to solve for the stream function by means 
of a series. 

In the first problem, the solution is carried far enough to afford an estimate for 
convergence. The torque on the inner cylinder and the second order steady streaming 
are computed, and a necessary re-examination of the meaning of the second concept 
is given. Various special cases of the second problem are shown to illuminate the con- 
cepts of virtual mass and viscous damping. The results for high frequency vibration 
bear on an experimental paper by Stuart and Woodgate, some of the low frequency 
results on theoretical papers by Andres and Ingard, and by Stokes. Both problems 
are of general interest in connection with the design of certain control mechanisms, 
and with experiments on Helium II. 

1. Introduction. Successful applications of conformal mapping to problems of 
fluid mechanics go back over a hundred years [3] and continue to the present. (See 
for example the many papers of L. C. Woods in the Proceedings of the Royal Society (A) 
and elsewhere.) Almost all of these applications have concerned two-dimensional flows 
of an inviscid incompressible fluid, since the well-known result 


V(x, y) , PU(x,y) _ EE n) , ve, 2] 
dx” oy ae on 


2 


dk | 
dz (1.1) 








V[x(E, n), yE, »)] = VE, n) 


concerning the transformation of the Laplace operator under a conformal transformation 
=29) @=x+ty,f =& + tm) 


means that the governing (Laplace) equation for the stream function remains invariant 
if the plane of the independent variables is subjected to a conformal transformation. 
Consequently, any awkward boundary shape can be changed into a more convenient 
one, the only penalty being (generally) an alteration in the conditions which the stream 
function must satisfy on that boundary. 

The object of this paper is to show that the use of conformal mapping provides 
an effective way to solve certain two-dimensional problems involving the flow of a 


~ 
« 
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viscous incompressible fluid. The governing equation for ¥ is now more complicated— 
(2.8) with J* = 1—and is no longer invariant under a conformal transformation, but 
if the awkward boundary is near a convenient one then the use of conformal mapping 
to exchange the former boundary for the latter will cause only a small alteration in 
the equation and boundary conditions, so that the transformed problem can be handled 
by standard perturbation methods. The particular problems discussed here involve 
flows forced by small amplitude motions of a boundary and therefore demonstrate 
the effectiveness of this approach in what might be termed “time dependent boundary- 
perturbation problems.” 

That conformal mapping would provide a useful method of attacking boundary 
perturbation problems was suggested several years ago by Lin [5]. In another paper [10], 
Segel has proved that the method is a valid one for such problems involving the equa- 
tion d°w/d2° + d°w/dy’ + Aw = O (dX a constant), e.g., for determining the possible 
modes and frequencies of vibrating membranes bounded by an ellipse of small eccen- 
tricity (“near-circle’’) or by a ‘‘near’’ isosceles right triangle. The present paper and [10] 
therefore comprise a detailed investigation, practical and theoretical, of Lin’s suggestion. 

The particular flow problems selected as examples are sufficiently complicated both 
to furnish a good indication of the scope and limitations of the method and to 
provide illustration and clarification of such general concepts of fluid mechanics as 
virtual mass and steady streaming. Furthermore, two-dimensional flows between moving 
nearly concentric circular cylinders currently are of great interest in connection with 
the design of control mechanisms for aircraft and rockets, and with experiments on 
Helium IT [14]. 

2. Rotating cylinders: formulation. Consider two right circular cylinders of radii 
r, and r, , r2 > 7, . Supposing that their respective axes are always parallel, we investi- 
gate the two-dimensional motion caused when the axis of the inner cylinder rotates 
in a circular path about the axis of the outer. Let the complex z-plane represent a cross- 
section perpendicular to the axes. The intersection of the plane with the cylinders is 
two circles; let z; = r,¢ exp [¢w a(¢)] and z, = O be the position of their centers at time 
t. Here 1/w is some typical time associated with the motion, ¢ is a dimensionless constant, 
and a is any real valued function of the time ¢. 

As is well known, the bilinear mapping ¢ = (Az + B)/(Cz + D) sends “circles”’ 
into ‘‘circles”’. (A “circle” is a circle or a straight line.) It would seem likely that proper 
choice of the constants A, B, C, and D would enable one to map the given eccentric 
z-plane configuration into concentric circles in the ¢-plane. This is indeed possible: 
if ¢ = p exp (i), the mapping 


. 2£+7,(6 — €) exp (twa) (2.1) 
S62 + (1 — 6e)r, exp (twa) ’ er 
with 
6 = —2f{(r3/i) —-1 —e€ + [{@e/r) -—1—-— &}? — 4€2]'74 7°, (2.2) 
sends the moving circle | z — z, | = r, and the fixed circle | z | = r, into the fixed con- 
centric circles p = 1 and p = 8. (Compare [4].) Here 
B = (r2/r,)(1 — € 6)". (2.3) 


Instead of employing r,/r,; and ¢, the actual dimensionless geometric parameters 
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of the problem, we will express the solution in terms of 8 and 6. To this end we note 
that 


e= —(8° — 1) 6/(1 — & 6). (2.4) 
In order to obtain the equations of motion in the ¢ plane, we observe that if u and v 
are the velocities in the direction of p increasing and of ¢ increasing respectively, then 
the continuity equation is [2, p. 114] 
O(pJ' *u) Op + a(J' *y) /db = (, (2.5) 
where the Jacobian J and the dimensionless Jacobian J* are given by 
| 2 22 
’ Z x, 4 ri—- 3d) . 
ee ee ee oe: Soe) (2.6) 
| dg | O(é, 7) (1 — 2 édpcos¢+ p 6) 
If we therefore define a dimensionless stream function y by 
u =rwp 'J~'*(dy/de), v= —riwd'*(dyv/dp), (2.7) 


the continuity equation is automatically satisfied. The equation of motion is now easily 
found. Writing the z-plane Navier-Stokes equation in vorticity form (i.e. eliminating 
the pressure) and using (1.1) and the related result 


OW, , V2.) ACY, , Wz) OE, 0) 


A(x, y) OE, n) A(x, y)’ 


we arrive immediately at the equation of motion 


d(Ay) 1 aA(Ay/J*, ) ] 
= — A(Ad/J*), 2.8 
ot p O(p, >) R A(ay¥ (2.8) 
where Ay = (d°~/dp) + p ‘(ay dp) + p (d°y dd’), R is the Reynolds number 
R = rw’y, (2.9) 


and v is the kinematic viscosity. 

Before establishing boundary conditions, we must specify the physical problem 
under consideration more precisely. For as the inner cylinder rotates as a whole about 
the origin, it may retain its orientation with respect to the line joining 2, and z, , or 
with respect to the outer cirele, or it may rotate about its own center in some other 
prescribed manner. We consider here the case in which the inner circle retains its orienta- 
tion with respect to the fixed outer circle. The other constant rotation cases require 
only a simple modification of the boundary conditions. 

The case in which both z, and z, are fixed, z, # 2. , and | 2, — 2, | is small, and in 
which each circle rotates about its own center with angular velocities w, and w, , is 
treated by Wood [15] who uses the time independent version of (2.1) appropriate to 
his problem, writing it down as a change of variables with no comment. He then expands 
the stream function and boundary conditions in powers of 6. Essentially the same 
problem as Wood’s was also treated by Nikitin [8]. The latter used polar coordinates 
and therefore had to expand the position of the boundary, as well as the stream function 
and the boundary conditions, in powers of the appropriate small parameter. 

In order to formulate our boundary conditions, we suppose, with no loss of generality, 
that a(0) = 0. Then if at ¢ = 0 a point on the inner circle has position 2(0) = rye + 
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r, exp (76), at time ¢ the point will be at z(¢) = rie exp (tw a(t) ) + 17, exp (20) and will 


have velocity z'(t) = ir,ew a’(t) exp (iw a(t) ) with components a’wr,e sin (6 — wa) 
and a’wr,e cos (@ — wa) respectively normal and tangent to the inner circle. But, from 
(2.1), when z = r,e exp (twa) + 1, exp (26) , 


exp (i¢) = {exp [i(@ — wa)] + 6}/{6 exp [7(@ — wa)] + 1} 


ft 
With this we can easily find cos (@ — wa) and sin (@ — wa) in terms of ¢: 
cos¢ — 26 + & cos¢ ’ : sin ¢(1 — 5) / 
Ah ele EB a sin (0 — wa) = Ss (2.10) 
1—2écos@¢@+6 ” ae ae 1—2é6cos¢+ 6 


The boundary conditions on y are therefore 


cos (6 — wa) = 











ay a’(B’ — 1) &1 — &) sing ; 
tp=1: —~= —- 7 rar wer X (2.11) 
mad ag (1 — B81 — 25 cos + 8)?’ 

dy _ a’(s* — 1) (1 — &)(cosg — 2 6 + & cosg) 

dp (1—p6 &)\(1 —2écos¢+ 6) 
at p = B: Oy/dd = OY/dp = 0. (2.12) 


As will be examined in more detail later, the boundary conditions (2.11) and (2.12) 
must be supplemented by the condition that the pressure be single-valued. Then the 
parabolic differential equation (2.8), the boundary conditions (2.11) and (2.12), and 
an initial condition of the form 


t = 0: ¥(p,¢, 0) prescribed, (2.13) 


together comprise a well-defined mathematical problem. 
The problem is non-linear. In order to be able to solve it, we restrict ourselves to 
between the centers is small: that is, to the 


the case in which the distance | z, — 2, | 
case in which e, and hence 4, is small. In the standard way, we assume that, for suffi- 
ciently small 6, y is given by the series 

¥= 7,67 2 Oo T °° (2.14) 
and we expand the differential equation and boundary conditions in powers of 6, set 
the coefficients of 6, 5°, --- equal to zero, and obtain 


0(6): AAy, — R d(Ay,)/at = 0. 





p = 1:4y,/dp = a’(B’ — 1) cos¢, dy,/d¢ = —a’(s’ — 1) sing. 
p = B: 0y,/dp = Oy, /dd = 0. (2.15a, b, ce) 
0(8’): AAy, — R d(Ay,)/dt = = Ayr » Av) + 4 A(p cos¢ Ay,). 
p Op, >) 
p = 1: 0y./dp = 2(8’ — 1) cos 24, Oy./d@ = —2(8° — 1) sin 2¢. 
p = B: 0y./dp = dy./dd = 0. (2.16a, b, c) 
| Ay.) RI av, , Ave) , Ave, Avs) , Avi, — 4p cosd Ae | 
0(5°): AAy,; — R——— = = | —— + 
Vs dt p a(p, >) a(p, d) ba a(p, >) 


+ A[4p cos @ Av, — 2(1 + p’ + 2p’ cos’ ¢) Ay, ] 


1: dy;/dp = (B" — 1)(8? — 4 + 6 cos 24), 
dy;/d¢ = —(B’ — 1)(8? + 2 + 6 cos 24). 


8B: dy;/dp = dY;,/dd = 0. ete. (2.17a, b, ce) 


> 
Il 


> 
Il 
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Once y is found, the velocities u and v follow at once from (2.7), and the pressure 
p(p, ¢) = bp, + Sp. + --- can be found from the Navier-Stokes equations in (p, ¢) 
coordinates (cf. [2, p. 101-103]). Since we will need it later, we note that p, satisfies 


po (Op,/Ap) = vwp *[A(A¥,) 4¢)], po (8p, /d@) = —vwp[A(Ay,)/dp]. —(2.18a, b) 


A formula for an important physical quantity is that for the counterclockwise torque 
about the origin which it is necessary to exert in order to keep the inner cylinder rotating 
in the prescribed manner. To derive it, we consider the moment caused by the stresses 
P,, and P,, acting on a given element Ag of the inner circle. With the aid of a simple 
diagram and a little trigonometry, the counterclockwise moment exerted by the fluid 


on the element is seen to be 
[Pre sin (@ — wa) + P,gri(1 + € cos (6 — wa)]J'”* Ad. 


The factor J'’* must be included since P,, and P,, are actual stresses; to find the force 


they exert on p = 1 one must therefore multiply by the length element | dz | = 
| dz/d¢||d¢| = J'” dd. To obtain the desired torque 7, we change the’sign and integrate 
7 —r; | {P,,esin (@ — wa) + P,,[1 + € cos (@ — wa))]}J'” do. (2.19) 


P,, and P,, can be written in terms of u and v, and hence y, by [2, p. 103] so that the 
coefficients 7, in the expansion T = 7, 6 + T, & + --- can be found using (2.10) 


and (2.4). 


3. Constant angular velocity. We consider now the special (and most interesting) 
case when the inner circle rotates about the origin with constant angular velocity w. 
Here a(t) = #, a(t) = 1. The boundary conditions are time independent, so that y 
is as well. The equation for y, is then AAy, = 0 so, from the boundary conditions, 
the form of y, is evidently 

v, = (8B — 1) cos ¢e(p). (3.1) 


Since e(p) satisfies an equi-dimensional ordinary differential equation, it is easily found 


to be 


e(p) = Cip ' + Cop In p + Cap + Cyp'. (3.2) 
From the boundary conditions (2.15 b, c), with a’ = 1 of course, we find 
B° ; — (8° + 1) 238° Ing+2Ing— 6s +1 l 
C,=- C, = “+4 1 ¢, = ee ==>, 
' s4,* ~* d, — 2d, aa $4) 
where 
d, = 6’ Ing+ing — 6B +1. (3.4) 


But the solution just found is not the only one satisfying AAy, = 0 and the boundary 
conditions (2.15 b,c). In facet 


h(p) = b, + b, In p + bsp’ + byp’ In p (3.5) 


satisfies AAh = 0, so provided only that the two requirements h’(1) = 0 and h’(8) = 0 
are satisfied, (8° — 1)(e)(cos ¢) + h satisfies the same equation and boundary con- 
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Fig. 1. Shapes of e(p)/p and e’(p). 
ditions as (8° — 1)(cos ¢)(e) itself: the four arbitrary constants in h are not determined 
by the two above requirements on h’. Recall however that the pressure must be single- 
valued. Since 0(Ah)/dp = 4b,p°', we see from (2.18b) that unless b, = 0, p, will have 
a term proportional to ¢ and will not be single-valued. Hence b, = 0, the requirements 
on h’ imply that b, = b; = 0, and the stream function is determined up to an arbitrary 


additive constant. 
The solution obtained is the only time independent one to the O(6) equations (2.15) 


with a’ = 1. However, any time dependent solution to (2.15a), having vanishing space 
derivatives on p = 1 and p = 8, could still be added to the y, of (3.1), and appropriate 
initial conditions could be satisfied by superposition of these solutions. The result would 
represent a flow forced by the assigned initial conditions and by the revolving cylinder. 
We are interested only in the motion caused by the revolving cylinder; consequently 
the y, of (3.1) is, as it stands, the desired first order stream function. 

Before considering the higher approximations, let us examine some of the features 
of the flow as given by the first order solution. We have performed calculations for 


8 = 2.718 corresponding to a configuration in which the outer radius is about three 
times the inner. The particular value of 8 was of course chosen in order that In 8 = 1, 
thereby simplifying the calculations. 

With 8 = 2.718, the functions e(p)/p and e’(p) have the shapes given in Fig. 1. 
Since the first order approximations to the velocities are u, 6 = (p dy,/d¢) 6 = 


(— 6)(8° — 1) p ' esin 6 and v, 6 = (— dy,/dp) 6 = (— 4)(6° — 1) e’ cos 9, the “flow” 
in the ¢ plane forms two symmetric vortices. At a given instant, the flow in the physical 
plane will therefore have the slightly asymmetric appearance of Fig. 2. 
Since the inverse of (2.1) is 
(1 — €6) te— 6 ita rae ; 
=f, exp (tw!) ——— ———— == F(t) exp (iwt) (3.6) 


l — of 
the velocity field looks steady to an observer rotating with angular velocity w. (This 
also follows at once from the fact that an observer sitting on the inner cylinder would 
see a fixed geometrical configuration in which moves a fluid driven by the steady counter- 
clockwise rotation of the outer cylinder.) The fluid is pushed away from the front of 














1961) VISCOUS FLOW BETWEEN MOVING CIRCULAR CYLINDERS 341 





‘ -->—-~ \ \ 
’ Sinaia” : 






--- ¢-—~ 
‘ 
‘~-p-—- 









~---g@- ---~ 






~-—-»--"~ 
Fic. 2. Rotating inner cylinder: instantaneous streamlines. 


the inner cylinder and then is deflected outward and backward by the outer cylinder. 
It flows inward and forward again at the back, filling the gap left by the retreating 
inner cylinder. 

We must explain briefly why we have spoken of “flows” in the ¢ plane. If one uses 
conformal mapping for determining the flow of inviscid incompressible fluids, it is 
certainly permissible to speak of a flow in the transformed plane, since the governing 
equation (Laplace’s) is invariant. Here, the governing equation is changed so that the 
“flow” in the transformed plane is not generally speaking a possible flow of a real fluid. 

If 5 is small enough to insure convergence of the y series, the flow should not depart 
appreciably from the first order flow described above. Nonetheless, it is necessary to 
calculate higher approximations in order to find how small 6 need be—and to compute 
such quantities as the torque 7’, for (3.1) can be shown to give 7, = 0. (This can also 
be deduced from symmetry arguments which demonstrate that 7’ must be an even 
function of 6.) 

Therefore we note that with the value of ¥, already computed, the 0(6°) equation 


(2.16a) becomes 


, - ‘ CC; — 8C,C 7 
AAy, = R(g’ — 1)’ sin 244 c2(2 — tt) + SaaS Ses 5 sci | 
p p 


+ (87 — 1)(64C, — 16C.p~” cos 2¢), 


(3.7) 


so that its solution is of the form 
y,. = (8? — 1)[f(p) cos 26 + g(p) sin 2 + h(p)]. (3.8) 


Substitution of (3.8) into (3.7) and solution of the resulting differential equations subject 
to the boundary conditions (2.16 b, c) show that 


f(e) = A,p~ + A, + Asp" + Ayp* + C.p° In p, (3.9) 
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g(p) = Bip’ + B, + B;p° + Byp* 


2 Cie? , OCs 
+ RG — of 48 g ? me 


Y 


tr 


16C,C, — C? — 4€.C; » Cz 2) 2 
mt :) 8 3 p In p— 16 p in | ; (3.10) 


(3.11) 








h(p) = D, + Dz, In p + Dsp’ + Dap In p+ a. 
By an argument parallel to the one around (3.5), it can be shown that p, is single valued 
if and only if D, = 0. Computation of the correct values for the remaining constants 
is now straightforward but tedious. The results are given in [9] as are the complicated ys 
equation and its lengthy solution. With them one can derive the torque formula 


T>/vpyer, = (6? — 1)[2C, + C2 + 2€,(48° — 1) — (8? — 14+ C, + 3C, — C,)). 


Also, since the convergence radius of the assumed expansions in powers of 6 cannot 
be estimated rigorously, the magnitudes of the higher order corrections to Y must be 
used to give a “‘practical’’ estimate. All we know a priori is that we have neglected the 
nonlinear convection terms, as opposed to the diffusion terms, in (2.8) so that we must 
have 6° “much less” than 6/R. Although we can choose R to be as large as we like and 
then choose 6 to make 6 < 1/R, we have assumed that diffusion is the dominant process 
in the flow so that our approximation is essentially a low Reynold’s number one. 

Calculations of the tangential velocity v [9] for 8 = 2.718 and R = 1 show that 
if 6 = — 0.01, the third order correction to v is at most 2 1/2 per cent of the second 
approximation so that in this case ‘convergence for | 6 | much less than 1/R” means 
“convergence for | 6 | less than about 10°*/R’’. If one used for a fluid “Rotring” motor 
oil (v = 9 em’ sec’, 1» = 8 gm cm” sec” at 1 atmosphere and 20°C according to p. 230 
of Schlichting’s Grenzschicht-Theorie) the case 8 = 2.718, R = 1, 6 = — 0.01 might 
correspond to an inner cylinder of radius 1.2 em rotating once per second within an 
outer cylinder of radius 3.2 cm, the axis of the inner cylinder being displaced 0.08 cm 
from the axis of the outer. This displacement is 6.4 per cent of the inner radius and 
4 per cent of the outer. The value of the torque per unit length in this case is 6.7 XK 10° 
dyne-cm/cm. 

Many details of the flow are most easily determined by obtaining the required 
quantity from ¥(p, ¢) evaluated at the point (p, ¢) corresponding to the point (r, @) 
in question: the torque calculation was carried out in this manner. In some instances, 
however, it is necessary to find W(r, @) explicitly by developing from an expansion 
of (2.1) series such as 

p=s-+ 4(8 — 8°) cos(@ — wt) +--- (s = r/r,), 
cos @ = cos (6 — wt) + 6(1/2)s"'(6° + s°)[1 — cos 2(0 — wt)] + ---, 


and thereby computing V(r, 0) = V(r, 0) + 6, (r, 6) + --- from Y(p, ¢). 

The most important quantities to be found from W(r, @)—or, more correctly, from 
V(r, 6, t)—are the steady streaming velocities u,(r, 6) and v,(r, 6). Since a computation 
shows that 

T 2 2 2 
v,=T" [ V(r, 6, i) dt = fe ° — e(s) + Ee = : e’(s) + no | + 0(6°) 


2 





= §V,, + 0(8) (T = 2r/w), 
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the leading contributions to the streaming are u, = 0 + O(6°) andv,/(r,w) = —dV2,/dr = 
—(C, + 4C,)(6°s"* — 8s) giving (to this order) a steady tangential flow opposite to the 
direction of rotation, its magnitude increasing from zero at the outer cylinder (s = 8) 
to — (C, + 4C,)(6° — 1) at the inner cylinder (s = 1). The fact that the direction of 
the streaming is opposite to what one would expect seems to be common to low Reynold’s 
number problems. (See Sec. 6.) 

It is at first alarming that although the boundary conditions require that the tan- 
gential velocity at a point (r, 6) of the inner cylinder be proportional to cos (@ — wf) 
and therefore have a zero time average, our calculation shows a non-zero streaming 
at the inner boundary. To resolve this difficulty, it is necessary to examine more closely 
the meaning of these calculations. What we have done is to compute the average velocity 
at a given point (r, @) from a coordinate system fixed in the stationary outer cylinder; 
the results could be verified experimentally by suspending an appropriate device at 
(r, 0) from a frame fixed to the outer cylinder. Viewed either mathematically or experi- 
mentally, it is clear that this procedure can have no meaning for points inside r = 
r,(1 + «), for every such point is for a time inside the rotating inner cylinder. During 
this time, the velocity at such points is undefined so a formal mathematical average 
has no meaning, and obviously an instrument, even if it could pass through the cylinder, 
would produce irrelevant readings. 

If one were standing on the inner cylinder, no difficulties in measuring the streaming 
would be apparent near it, but the same difficulties would appear near the outer. The 
theoretical result here is u, = 0, v, = riw(C. + 4C,)(s — 1), again giving a tangential 
streaming opposite to the direction of rotation seen by the observer. To sum up, in a 
flow caused by boundaries in motion relative to each other, there is a different steady 
streaming with respect to each of the boundaries; to single out this aspect of the flow 
generally has no meaning near to boundaries other than that one selected for reference. 

4. Variable angular velocity. If the inner cylinder rotates with a variable angular 
velocity, i.e. if a is no longer considered to be a linear function of ¢, we write the first 
order equation (2.15) as 

0Q,/at = R™ AQ, , Q, = Ay, . (4.1) 
If we assume separable solutions of the form 
Q, m(p) (8° — 1) cos @ exp (ipi), vy, = n(p)(8’ — 1) cos¢ exp (tipi), (4.2) 


then equations (4.1) are satisfied, for any constant p, if and only if 


n+ p'n’ — pn =m, (4.3) 
pm’ + pm’ + (k’p’ — 1)m = 0, k? = —ipR. (4.4) 

Since (4.4) is a Bessel equation 
m(p) = AJ,(kp) + BY, (kp); (4.5) 


n(p) can now be obtained by variation of parameters. In general, n(p) will contain two 
more arbitrary constants (in addition to A and B); in this case satisfaction of the two 
boundary conditions (2.15c) at p = 8 requires these constants to be zero. Since p is 
arbitrary, we superpose the solutions already obtained and attempt a solution of the 
form 


¥, = (8° — 1) cose [ exp (ipt)[oF(p, k) — p'G(p, k)] dp (4.6) 
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in which k = (— ipR)'” 
8 


F(p, k) = | [A(k) J (hE) + B(k) Y,(ké)] dé, (4.7) 


»8 
G(p, k) = | [A(k) J (ke) + Bk) Y (ke) Je dé, (4.8) 


or, since the integrals in (4.7) and (4.8) can be evaluated, where 
F(p, k) = —k "A(k) [J o(kB) = J (kp) | = k 'B(k) [ ¥ o(kB) —_ Y.(kp)], (4.9) 
G(p, k) = k' A(k)[B°J2(kB) — p’Jo(kp)] + k “BCA [8° ¥2(k8) — p’ Y2(kp)). (4.10) 


From the two remaining boundary conditions we therefore find that 


; kF(k) .., reaver ay _ _ KF(R) rr 
A(k) = Dik) [Y.(k) — B’Y.(6k)], B(k) = Dib [.J.(k) B’J2(Bk)], (4.11) 
in which 
D(k) = [Jo(k) — Jo(Bk)|[Y2(k) — BY .(Bk)] (4.12) 
— [J.(k) — BJ.(Bk)][Yo(k) — Y,(Bk)]. 
With all this, we see finally that y, can be written as 
a\~1/87 08 r° svt} PZi(p, k) — p 'Zo(p, k) |, 
y, = (2r) "(8 — 1) cos $ | _ | — — H(p) dp (4.13) 


where we have made the abbreviations 
Z(p, k) = [J.(k) — BJ (Bk) ][Yo(Bk) — Yo(pk)] 
— [J.(Bk) — Jo(pk)][Y.(k) — 8 Y.(Bk)], 
Z.(p, k) = [Y.(k) — B’Y.(Bk)][B°J2(Bk) — p J2(pk)] 
— [8° Y.(Bk) — p Y.(pk)][J.(k) — 8 J.(8k)], 


and where H is the inverse Fourier transform of a’: 


a’ (t) = (2r)7'”" | e '’’ H(p) dp. (4.14) 
(4.13) may be partially checked by noting that it should reduce to (3.1) when a’(¢) = 1. 
But when a’ = 1, (27)~'”H = 4(s), the Dirac delta function, and evaluation of (4.13) 


is afforded by finding the limit of Z,(p, k)/D(k) and Z.(p, k)/D(k) as k — 0°. This is 
most easily accomplished directly from the series definitions of the various Bessel func- 
tions involved; the computations show that (4.13) does reduce to (3.1) when a’ = 1. 

Finally, behavior of y, at large time can sometimes be found at once from known 
asymptotic formulas for Fourier transforms. Using a formula in ([7], p. 483), since for 
p — 0 the square bracket in the integrand of (4.13) behaves like 


1 (8? — p)(3p — p ') — (8° + 1) In (B/p) . 
— ; ~——., - O(p), (4.15) 
2 6B —1—(e +1) lInB + OG 


H (p)~Cp”" for p—O (where H,(p) = 0 for p <0) 
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then asi— ~ 
; CT(m + 1) (1/2)(8° — p°)(8p — p™') — (6* + 1) In (6/p) 
n~ es =f) ae] 2 2 ; . 
t f-i-Y+) ms 

The 0(p) term in (4.15) has been evaluated, but it seems too complicated to be of interest; 
note, however, that terms of 0(p In p), which arise from the Bessel functions of the second 
kind Y, and Y, , cancel, so that the remainder in (4.15) is O(p). 

5. Vibrating cylinders: formulation and first approximation. Let us again consider 
the two circular cylinders of Sec. 2, only this time let the center of the inner cylinder 
= 2, = 0 while z, , the center of the outer, vibrates about r,¢ according 


“£2 





remain fixed at z 
tO Z2 re (1 + cos wt). Neither cylinder will be allowed to rotate about its center, 
though, as before, such rotations present no difficulty aside from more complicated 
computations. This problem is a natural complement to that of Sec. 2; it also contains 
different features which may be associated with the fact that there is no longer a co- 
ordinate system in which the problem is steady. 
To solve, we proceed as before. The conformal mapping 
¢ = (e—7, d/(r, — 6&2) (5.1) 


transforms the circle with center z, and radius 7; into p = 1, and the circle with center 
z, and radius r, into the concentric circle p = 8. Note that (with the abbreviation a = 


r,/r,) the value of 4, 
2e(1 + coswt) 


sii ~ a? —1—€(1+ coswt)’+ | fa’ —1—e(1 + coswt)”}?—4e(1 + coswt)?}?/? (5.2) 





varies with time so that the mapping (5.1), and the Jacobian 
J = | dz/d¢ |? = ri(1 — &)°(1 + 2p 6 cos + 8p)” (5.3) 


are now functions of time. The dimensionless stream function (defined as in (2.7)) 


will therefore satisfy 


A ay 1a(Ay/J, py) 1 ; a 
: oe SE oe J =w 5. 
} OT ( J ¥ p O(p, >) R A(ay/J) (7 wed 6 4) 
with appropriate boundary conditions on p = 1 and p = 8. The outer circle, while 
now concentric with the inner, has a time varying radius 
8 = all — de(1 + cos 7)]’. (5.5) 


To make the boundaries into fixed concentric circles, one would have either to make 
a further radial (non-conformal) transformation or to use the standard series expansion 
method, discussed in [10]. However, as 


B=ate(1 + cos 7)’alo’® — 1)' + --- = afl + O})], (5.6) 


the effect of the variation of 8 with*time need not be taken into account until the third 
approximation is considered. 
Taking a and ¢ as the geometric parameters of the problem, we assume that y = 
y, + ye’ + --- and proceed as usual. Noting that 
—j = e(] a. cos t)(a’ _ no + 2*(1 + cos t)*(a’ = 1)? + ital (5.7) 
we find that (5.4) separates into 


O(e): R A(Ay,)/dr — AAYy, = 0. (5.8) 
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4k p cos ¢ sin t AY, 
* (5.9) 
+ Al4p cos ¢ (1 + cos 7) Ay, ](a’ — 1)", 


O(e’):R d(Ay2) ‘or = AAy, = —Rp" (Ay, , ¥:)/A(p, @) + 


while the boundary conditions are 
O(e) : p = 1: dy, /dp = Ay,/dd = O. 
p = a: 0y,/dp = —sing@sin 7, dy,/dd = —a cos@sin r. (5.10) 
O(e) : p = 1: dy2/dp = Ay2/d¢ = 0. 


p = a: dy./dp = —ala’ — 1)7'(2sin r + sin 27) sin 2. 


dy./d6 = —a’(a’? — 1)7'(2sin r + sin 27) cos 29. (5.11) 
As (5.8) is the same as (4.1), we follow the procedure of Sec. 4 to obtain the following 
solution to (5.8) and (5.10): 
¥, = Im [e(p) sin ¢ exp (77)], (5.12) 
where, if S = 7°?R’” 
e(p) = —(A/2S){p[Jo(Sp) — Jo(S)] + p'[eJ2(Sp) — J2(S)]} 
— (B/28){p[Yo(Sp) — Yo(S)] + p '[e Y2(Sp) — Y2(S)}}, (5.13) 
A = 2S[a’Y,(Sa) — Y.(S)]/d(@), B = —2S[a’J.(Sa) — J(S)]/d(a), (5.14) 
d(a) = [J.(Sa) — Jo(S)]la’ ¥2(Sa) — Y2(S)] 
— [a°J,(Sa) — J2(S)][¥o(Sa) — Y,(S)]. (5.15) 
The most interesting physical quantity to be derived from this result is Fy , the 
horizontal force on the inner cylinder caused by the prescribed motion of the outer. 
(The vertical force is zero by symmetry.) A computation shows that the force (taken 
positive in the direction of x increasing) on a length L of the cylinder is 
Fy, = Im {R7'MUw exp (iwt)S[AJ.(S) + BY.(S)][1 + O(6]}, (5.16) 
where, for reasons discussed below, we have made the abbreviations 
Uso = nr, €, M = porriL. (5.17) 
The problem (call it P) under consideration is closely related to the problem P’ 
of determining the motion caused if the inner cylinder vibrates (z, = — r,e(1 + cos wt)) 
while the outer remains stationary (z. = 0). The connection between the two is seen 
by the following standard argument. The problem P originally discussed is governed by 
du/dt + u du/dx +0 du/dy = —po' dp/dx + r(Gu/dx” + d’u/dy’), (5.18) 
the corresponding equation for v, the continuity equation, and the boundary conditions 
appropriate to the vibration of the outer cylinder. (For this paragraph only, u and v 
represent velocities in the x and y directions.) If we make the change of variables 


, 


z’ = 2 — er,(1 + cosa), ’ = tft, y’ =y, uw =utnr,ewsinal, (5.19) 


and 
p’ =p — pew coswt X= p— por, ew coswt r cos 8, (5.20) 
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we get a new problem with equations identical with (5.18) etc. but, due to the change 
from x to x’, with boundary conditions appropriate to problem P’. Consequently a 
solution to P’ can be related to a solution of P by (5.19) and (5.20). In particular, the 
pressure on the inner cylinder for P’ is found from the corresponding pressure for P 
by subtracting por,ew cos wl cos 6. This means that the force on the inner cylinder 


in P’ is found from (5.16) by adding 
A’ = MU w cos wt. (5.21) 

Until now, the only assumption we have made is that | ¢ | < 1, i.e. that the amplitude 
of vibration is small compared with the radius of the inner cylinder. As is seen from 
(5.4), this allows us to neglect convection terms (to a first approximation) regardless 
of the magnitude of R. That is (5.12) and (5.16) give the first order solution for small 
amplitude vibrations, regardless of their frequency. As these expressions are complicated, 
it is however advisable to examine their behavior for high and low frequencies. 

or large R (and hence large | S |), ¥, = Imfe.(p) sin @ exp (ir)] where, using the 
asymptotic formulas for J and Y 


9\'/2 19 oe 1 T 
~ Ss K 5 7 *) r ( ~ 1) ( = *) 
€..(p A. (2) 2 p sin (s 4 p ; cos | S 4 
- a (, + 18) sin (s — *) 
kee t ae RS T 
— , rS S : ; (s - *) ( ™ 1) ; (s ~~ *) 
I (2) 2 p cos | Sp — 3 +\p ’ sin|S r 


_ a (, 4. 3) cos (s - *) ‘ (5.22) 


A. = dz'(2/r8)'"[sin (S — 2/4) — a®” sin (Sa — 2/4)], 





B.. = dz"(2/r8)'”"[cos (S — 2/4) — a®” cos (Sa — /4)], (5.23) 
d. = (2/rS)a~'"(a* — 1) sin [S(@ — 1)]. 


It is necessary to include the “‘small’’ O(1/S) terms in the square brackets of (5.22 
because the “‘order one’’ term is in fact zero at the inner boundary. The region in which 
the O(1/| S |) terms are important is given approximately by (p — p') < | 1/S | 
orp <1+R'?,r = rnp < 1+ (v/w)'”. As the term proportional to (p — p™') is 
the inviscid solution, we see that for high frequency oscillations the influence of vis- 
cosity is primarily confined to a boundary layer near the inner cylinder having thickness 
of order (v a)” 

Using (5.22), 
due to high frequency oscillations of the outer. The result is 


one can compute the asymptotic form of the force on the inner cylinder 


Fw = a'(a® — 1)7'MUww{ [2 + 2(2/R)'”] cos wt + 2(2/R)'” sin wt}. (5.24) 

From (5.21), when the inner cylinder is oscillating with high frequency as described 

in Problem P’, the force exerted by the cylinder (i.e. the force required to move it) 
is then 

—Fhe = a'(e? — 1)7'MUww{[a*(a? + 1) + 2(2/R)'”] cos wt + 2(2/R)'” sin wt}. (5.25) 
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Since the center of the inner cylinder is at r,«(1 + cos wf), it has velocity V = r,ew sin 
wt = U, sin wt. Consequently, if a length L of the inner cylinder had mass M* it would 
require a force of M/*wU> cos wt to move it in the absence of a fluid. The fluid’s presence 
causes the appearance of an additional force MU w(a*” — 1)~* [a*® + 1 + 2a°(2/R)'””] 
in phase with the acceleration, so the cylinder behaves as if it had the additional “virtual 
mass” M(a* — 1)" [a® + 1 + 2a°(2/R)'*]. The dependence of this expression on 
a reflects the influence of the fixed outer cylinder on the force required to move the 
inner. If the outer cylinder is made to recede indefinitely (a — ) the virtual mass 
approaches M[1 + 2(2/R)'*], which is Stokes’ result for a vibrating circular cylinder 
in an infinite fluid [12]. As R — © the virtual mass approaches M, the mass of fluid 
displaced; as it should be, this is the inviscid result. 

The concept of virtual mass arises from the fact that a body moving in an inviscid 
fluid is subject to an extra inertia force which can be found by augmenting the body’s 
mass by a multiple of the mass of fluid it displaces, the numerical value of the multiple 
depending on the nature of the motion and the shape of the body. If the fluid is slightly 
viscous, the inviscid solution still holds, approximately, except for the boundary layer 
region of slow motion near the body, so that one expects to find a correction to the 
virtual mass due to the increase in the apparent size of the body caused by the boundary 
layer. (We assume that separation does not occur.) In the problem considered here, 
the viscous correction term to the virtual mass, [2a°(2M)'*]/[(a° — 1)R'*], may be 
interpreted in this way by writing it as 2'’’Ap.La*/(a” — 1) since A = 2z7,(v/w)'” 
is of the order of the area effectively added to the inner cylinder by the boundary layer. 
In addition to the virtual mass force, there is a viscous damping force 


—2(2 R)' “a (a® — }) "MU ,w sin wt 


on the body, opposing its motion. The coefficient here approaches 2-2" 2.83 as 
a — ©, In experiments by Martin [6] and by Stuart and Woodgate [13], the values 
2.87 and 2.94 were found for this constant. The small discrepancies could be due to 
any of several factors, but as some finite container for the fluid surrounding the vibrating 
cylinder must have been used, it is worth noting that a bounding outer cylinder of 
8.5 and 5 times the radius of the inner could have produced the observed discrepancies. 
Stuart in fact recalls that the ratio a for their experiment was of the order of 10 (the 
outer container was cylindrical) so that about one-fourth of the observed discrepancy 
is attributable to the presence of the outer cylinder. Martin used a vibrating wire in 
a complicated set up: a photograph of his apparatus (in his paper) indicates that the 
nearest obstacle is so many diameters away from the wire that the effect discussed 
here is insignificant in his experiment. 
For small R, the first order stream function is approximately 


¥, = (2d,)~™' sin g@sin wi[p + (a —1)p—a'p ' — Qa? + 1)pln pl, (5.26) 
7) 


‘a4 


i, = (a? + 1)Ina— a’ + 1, (oO. 


bo 


so that the response of the fluid is very nearly in phase with the forcing movement of 
the boundary. The leading term in expression (5.16) for the force on the inner cylinder 
due to vibration of the outer is 


Fo = — 4(a’ + 1)(d)R)"' MU ww sin wt; ‘5.28 
(9.28) 
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since A’ of (5.21) is independent of R, this is also the leading term in the expression 
for the force on the inner cylinder in problem P’. The principal force in low Reynolds 
number vibration is therefore a viscous damping which becomes infinite as R approaches 
zero. 

Since d[(a® + 1)/d.]/da = — (a® — 1)’a' dj? < 0, the viscous damping decreases 
from an infinite value as the ratio of the radii increases from unity (at fixed R). Note 
however that (5.28) is valid only for finite a as it was derived under the assumptions 
that | S| and | aS | are small. That is, mathematically, a small value of | aS | was assumed 
in obtaining simplified expressions for A and B in (5.16) so that in discussing the variation 
of (5.28) with a, we assume that RF is fixed, and so small that (5.28) is valid for a range 
of a. The solution for small Reynolds number oscillations in an unbounded fluid can 
only be obtained by first letting a tend to infinity and then taking the limit for small R. 
(See Stokes’ results described below.) Physically, if | S| = R'” = r,/(v/w)'” is small, 
then the thickness of the boundary layer is large compared to the radius of the inner 
eylinder: if | aS | = r,/(v/w)'”* is also small, then the thickness of the boundary layer 
is also large compared to the radius of the outer cylinder. Consequently, in deriving 
(5.28) we have assumed that viscosity is important in the entire region between the 
two cylinders. 

Taking a further term in the expansion of (5.16) for small | aS |, we find that the 
first correction to (5.28) is a force 

‘ 5.4 4 9.2 » 2 
a + 3a ba In re Jar — ha In ow 5 MU.» cet = —Ci(o) MU, outed (5.29) 
acting on the inner cylinder due to the vibration of the outer. This means that for small R 


the virtual mass coefficient for problem P’ is very nearly C)(@) — 1. Examination of 
(5.29) shows that Cy(a@) is infinite fora = 1 and fora = ©, and therefore has a minimum 


for some intermediate a. 
As in [12], Stokes’ result for the force on a cylinder oscillating at small R in an un- 


bounded fluid can be found from the general formula of linear theory 


( r 17 -p\1/2 
Py = Im\- MU ge - Pag Rapes |p 
Im {—t¢MU we'*'[K — iK’]} = —MU[K coswt + K’ sin wt], (5.30) 
where K, is the modified Bessel function of the second kind as defined in [16], p. 374]. 
For small R, K and K’ reduce to 
K =1+7/R(L’ + x°/16), K’ = 4L/R(L*? + 2°/16), (5.31) 
L = In (4/R)'” — vy, y = Euler’s constant = 0.577 --- , 
which are the expressions obtained by Stokes [11]. If the oscillating cylinder is enclosed 
by another cylinder, our work gives the results 
K=Ca)-1, K’ = 4@’ + 1)/(@R). (5.32) 


A comparison of (5.31) and (5.32) shows that as R — 0 the damping coefficient becomes 
infinite like 1/(R In R) for a cylinder oscillating in an unbounded fluid, and like 1/R 
if the fluid is bounded externally, within the boundary layer, by another cylinder. 
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In the first case, the virtual mass coefficient also becomes infinite, more slowly, like 
1/(R In’ R) while in the second, the virtual mass coefficient remains finite. 

It is possible to find suggestive physical interpretations for the derived variation 
of the force with R and with a. Let us refer to the Stokes problem and the problem 
considered in this paper as the “unbounded” and ‘“‘bounded”’ problems respectively, 
and let us consider the force on the inner cylinder due to its own vibration. That for 
fixed a, in both cases, this force increases without limit as the Reynolds number tends 
to zero can be regarded simply as reflecting the fact that vibrations in a very viscous 
fluid are heavily damped. 

That as R — 0, the virtual mass tends to infinity in the unbounded problem but 
approaches a constant in the bounded problem can be illuminated as follows. In the 
solution to the unbounded problem, the radial factor in the stream function is a linear 
combination of 1/r and K,[r(zR)'’*] [12]. K, arises because of the inclusion of the dy/dt 
term and, mathematically, allows one to satisfy the viscous no-slip boundary conditions 
at the cylinder. Physically, as in [11], if Stokes’ solution is regarded as giving the behavior 
at large time of a fluid agitated by a circular cylinder which at some instant begins 
a sinusoidal oscillation from rest, then the dy/d¢ term can be seen to provide the inertia 
needed to curb what would otherwise be a continual increase in the amount of fluid 
carried by the cylinder as it persisted in its oscillations. (Similarly, if a fluid with constant 
velocity U at infinity flows past a cylinder, this curbing is accomplished by non-linear 
terms like u du/dx-linearly approximated by Oseen as U du/dx.) Since if R is very 
small, the velocities arising from K, are important until r becomes very large, it is 
clear that at small RF the influence of viscosity extends far into the fluid from the cylinder. 
The fact that the virtual mass in the unbounded problem tends to infinity as R — 0 
is thus seen to be due to the fact that the amount of fluid pushed and dragged by the 
cylinder—due to viscosity—also tends to infinity as R — 0. 

If the fluid is bounded externally, the situation is different. The radial coordinate 
is now bounded so that if F# is sufficiently small, the inertia terms are uniformly small 
compared to the terms retained. As we have seen, the leading contribution to the stream 
function is then the quasi-steady (5.26) which represents the in-phase response of an 
inertia-less fluid to the motion of its boundaries. The next most important term in the 
stream function gives rise to a finite virtual mass because the amount of fluid pushed 
and dragged is at most the finite amount contained between the two cylinders. 

From a simple point of view, the variation of the force with a should be governed 
by the fact that the nearer the outer cylinder is to the inner, the greater should be the 
influence of the outer cylinder, i.e. the greater should be the force on the inner—which 
accounts for the decrease of the viscous damping with increasing a (at fixed small R). 
But although the virtual mass first decreases as @ increases, it reaches a minimum 
and then starts increasing again. This is because the virtual mass is proportional to 
the mass times the acceleration of the fluid pushed and dragged by the inner cylinder. 
Now the closer the outer cylinder is to the inner, the faster u and v must change from 
the values imposed by the vibrating inner cylinder to the zero values imposed by the 
stationary outer cylinder, so that as a increases the derivatives of u and 2, i.e. the accelera- 
tions, decrease. But when a increases the amount of accelerated fluid increases. It appears 
that at small a the former process dominates, at large a, the latter. 

Once again, the above discussion of the qualitative effects of the outer cylinder is 


. “oe | 1] 1/2 2 1/2 s_ . . . : 
valid only if | aS | = aR’ = (ryw/v)'”* is small. Thus, discussing the pendulum experi- 
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ments of Bailey, Stokes [11, pp. 83-7] uses (5.31) in computing the correction to the 
virtual mass, and hence to the period, of the pendulum due to the thin wire by which 
the spherical pendulum was suspended. This is justified (and our calculations assuming 
small | aS | are not relevant) because although | S | ranges from 10~' to 10~? for the 
various cases considered, and is certainly small, the outer cylinder is so far from the 
wire that | aS | is of the order 10° and is thus not small but large. Appropriate simplifi- 
cations of (5.16) for large | aS | and small | S | (the case of a boundary layer thin compared 
with the radius of the outer cylinder but thick compared with the radius of the inner) 
could easily be worked out, but the effect of the outer cylinder on Stokes’ already small 
corrections would be insensible. 

6. Second approximation: small R. y, is found by solving the differential equation 
(5.9) under the boundary conditions (5.11); it clearly has the form 


Yo = sin 26/C,(p) + C.(p) sin r + C;(p) cos r + C,(p) sin 27 + C;(p) cos 27]. 


Since the solutions of the homogeneous equations are known, the inhomogeneous equa- 
tions for the C’s may be solved by variation of parameters. Finding y, therefore presents 
no difficulties other than computational ones, but these are formidable due to the com- 
plicated expression obtained on substituting y, into the right hand side of (5.9), and 
due to the fact that two of the four solutions to the homogeneous equations are also 
complicated combinations of Bessel functions. 

Interesting results are still forthcoming if we make the simplifying assumption 
that R is small and expand the known function y, and the unknown function y, in 


powers of R 


Vi = Vio +t Ryn + +++ , v2 = Yoo t+ Ry + ee. (6.1, 6.2) 
We obtain from (5.9) and (5.11) 
—AAYoo = (a* — 1) A[4p cos (1 + cos 7) Ayo], (6.3) 
p = 1: do0/dp = Ap2/dd = 0. 
9 = a: OWe0/dp = —ala’® — 1)7'(2sin r + sin 27) sin 24, (6.4) 


OWoo/0 = —a’(a’? — 1)7'(2sin r + sin 27) cos 29, 





a a 1 (Adio » Hoe) 4 Af4p cos eit cos t) Ay; ] 
p O(p, >) a 1 


Avro — Mabe, (6.5) 


4p cos ¢ sin r 
a —1 


a 


p = 1: dy2,/d8p = Ov2,/dd = O. Pp = a: On /Ip = Ayv2,/dd = O (6.6) 
from which we can obtain Yo and y2, . On substitution of Yio from (5.26), (6.3) becomes 
—AAYo = —8(a® + 1) do'(e? — 1)7"[2 sin r + sin 27]p’ sin 29. 

The solution of this equation and the boundary conditions (6.4) has the form 
Yoo = f(p) sin 26[2 sin r + sin 27], 
f(p) = A,p-? + Ay + Asp? + Acp* — (a? + 1)[2(e” — 1) do]~’p’ In p, 
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Fia. 3. Vibrating outer cylinder: steady streaming. (Since points in the hatched region are for a time 
outside the outer cylinder, a steady streaming velocity cannot be usefully defined in this region.) 


so that the equation for y¥., can be derived and the solution found to be of the form 


(cos 27 — 1)sin2¢, ,. 1+ 2cos7+ cos2r. 
vo = — a —h(p) - SS 851 26 h(p) 
d (a —1)d, 


+ 8(cos 7 + cos 27) sin 2¢ h;(p). 


The details are given in [9]. 
We observe the appearance of a second order steady streaming with stream function 


eR sin 26 l ; ; . =z 
ve. = 3 | a(o + ., + i(p — lef — (a& — 1)h, - i. 
Qa — ] p 


The first two terms arise from y, just as in the discussion of streaming at the end of 
Sec. 3. (We have expanded the e(p) in (5.12) as e(p) = eo(p) — iRe,(p) + O(R>).) The 
last two terms arise from the aperiodic part of ¥. (not present in the case of Sec. 3). 

A computation shows that for a = 2.718 the steady streaming forced by a slowly 
vibrating outer cylinder (Fig. 3) shows little qualitative difference from that forced 
by a stream slowly oscillating at infinity. (For the latter, see [l, Fig. 4].) One must 
beware of a (false) conception of steady streaming as some sort of actual flow; this 
would lead to the erroneous intuitive idea that the presence of the outer cylinder should 
cause the y,, streamlines to bend around into closed paths. In fact, as discussed in Sec. 3, 
the steady streaming velocity is an artificially singled-out part of the actual velocity. 
Failure to realize that this singling-out has no meaning near the outer cylinder (to be 
precise, in the hatched region of Fig. 3) would lead to the incorrect ideas just mentioned. 
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BOOK REVIEWS 
(Continued from p. 334) 


Wiener-Hopf idea, to treat the Summerfeld half-plane diffraction problem. The reviewer strongly 
supports the author’s view that here, and in other problems whose basic formulation includes a partial 
differential equation, the direct attack on the differential equation is far less laborious than the alter- 
native procedures which require the choosing of a Green’s function and the formulation of an integral 
equation as intermediate steps. This is especially true when the differential equation is not one of those 
few for which the appropriate Green’s function is known a priori. Also included in Chapter II is an 
outline of the treatment of several problems whose geometry and analysis closely parallel those of the 
diffraction problem. 

Chapter III is concerned with the details of the application of this technique to other problems 
in which the geometry is somewhat more complicated. Typical of these are the diffraction of waves 
by a pair of semi-infinite planes and the transmission of sound from a semi-infinite pipe. It is to be 
noted that the method becomes no more difficult in principle for such problems but the algebraic labor 
increases enormously. 

The remainder of the book is concerned primarily with variations on, limitations of, and approxi- 
mations within, this technique. The reader will find that he needs only Chapters I and II to understand 
both the basis of the method and the manner in which it can be used. The remaining material, including 
the rather complete bibliography, gives an excellent view of its exploitation and range of applicability. 
The book should be useful to those who need an introduction to the subject and to those who want 


a compendium of its accomplishments and limitations. 
GeEorGE F, Carrier 


Nichtlineare Mechanik. By Hans Kauderer. Springer-Verlag, Berlin, Géttingen, Heidel- 
berg, 1958. xi + 684 pp. $15. 51. 


In the book under review, the author has produced the type of volume which would be described 
in the English literature as a treatise on the vibrations of nonlinear systems. However, in addition, 
one aspect of nonlinear elasticity is treated as well. 

The first part (amounting to less than one quarter of the total volume), deals with elastostatics 
under nonlinear stress-strain relations but infinitesimal strains. A section dealing with theoretical aspects 
of elasticity theory under nonlinear (single-valued) stress-strain relations is followed by a discussion 
of problems of plane strain, torsion, and bending of thin plates. Nonlinearities arising from sources 
other than the stress-strain relation are not considered. 

Except for a short final section, the remainder of the book is devoted entirely to the vibrations of 
systems which are governed by ordinary nonlinear differential equations. More than two-thirds of this 
portion deals with systems having a single degree of freedom; the remainder contains a discussion of 
multi-degree-of-freedom systems. 

A feature of the book which strikes the reviewer as particularly pleasing is the logical arrangement 
of the subject matter. 

Vibrations in a single degree of freedom are subdivided into autonomous and non-autonomous 
systems. In the section on autonomous systems one finds, in that order, conservative, damped, and self- 
excited systems. The section on non-autonomous systems deals with the forced harmonic and sub- 
harmonic vibrations in the Duffing problem, with forced vibrations of systems capable of self-excitation, 
and with auto-parametrically excited systems. The stability theory of autonomous systems is an integral 
part of the section treating the response of these systems, while the stability of non-autonomous systems 
is treated in a special section. 

The last ninety pages are devoted to the vibrations of multi-degree-of-freedom systems. This 
portion is divided into systems having two, and having an infinity of degrees of freedom; each is again 
subdivided into autonomous and non-autonomous systems. The treatment of dual-mode non-autono- 
mous systems contains a discussion of forced vibrations of (otherwise) conservative, and of damped 
systems under simple harmonic excitation. The section of infinitely many degrees of freedom consists 


(Continued on p. 874) 
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A FUNDAMENTAL PROBLEM OF NAVIGATION IN FREE SPACE* 


BY 
GEORGE M. EWING 
University of Oklahoma Research Institute and U. S, Army Artillery and Missile School 


1. Introduction. We are concerned with the problem, termed fundamental by 
D. F. Lawden [6, p. 171], of transferring a rocket-driven vehicle from an initial position 
and velocity to a terminal position and velocity in such a way as to minimize the mass 
ratio. 

Mathematical models for this and other questions considered in recent publications 
(see Ref. [4] through [9] for examples) suggest problems of Bolza in the calculus varia- 
tions but often are not covered by the formulation of Bliss [1, p. 189] as a result of the 
essential nonincreasing character of mass. 

Treatments fall into: Type 1, in which mass M/ is permitted to have discontinuities 
(as in [4, 5, 6]) and Type 2, in which a bound is placed on the rate of mass flow (as in 
[7, 8, 9]). Type 2 simulates an obvious limitation of rocket motors. Type 1 problems 
though less realistic in this respect are of interest per se, sometimes appear easier to 
attack, and furnish bounds on the performance provided by solutions of corresponding 
Type 2 problems. Type 2 problems have been treated (see [9, p. 9, items (2), (3)]) with 
the aid of side equalities equivalent to the restriction that mass M be monotone, a 
device introduced by F. A. Valentine ({10, p. 4]) under more classical hypotheses. An 
alternative technique ([8, p. 582]) termed parametric representation is due to Miele. 
Such approaches bring Type 2 problems closer to the classical fold. In contrast it appears 
that definitive treatments of Type 1 problems require new lines of attack. 

The theory of optimal burning and steering programs is presently in an incomplete 
state. Many writers consider only the formal manipulation of first order conditions. 
There is no explicit statement of problems; hence the reader must decide from the 
context which functions are admissible and which steps are legitimate. Questions of 
existence and sufficiency are largely unmentioned. Miele ({[7, pp. 105-109] and [8}) 
gives sufficiency arguments for certain problems. 

In the formulation of Type 1 problems it seems desirable to admit general non- 
increasing mass functions M. Neither M nor the velocity is then necessarily AC (abso- 
lutely continuous) in the time t. Conclusions based upon differential equations or other- 
wise involving the integration of M’ or of accelerations can be correct only by coincidence. 

The present paper establishes the existence of the absolute minimum for a Type 1 
planar version of the problem mentioned at the beginning and shows that at least one 
optimal program is determined by certain solutions of an elementary system of equations, 
(3.2). The problem is partially covered by [4] and by the heuristic discussion of [6, 
pp. 171-175], subject to the objection of the preceding paragraph. 

2. Formulation of the problem. Let 7 denote the time required for a particle to 
go from point (0, 0) and velocity (u, v) to position (a, b) and velocity (U, V). Time T 
is preassigned except for Sec. 8. We understand 7’ to be positive with the single exception 
that the trivial case 7’ = 0 is included under Theorem 5.4 as type (v). 


* Received July 21, 1959; revised manuscript received November 12, 1959. 
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An admissible burning program is any real function M of time ¢ defined and non- 
increasing for all real ¢ and having respective constant values M(O—), M(7'+) for 
t<t > f. 

An admissible steering program is any real function @ of ¢ defined and continuous 
for all real ¢. Continuity appears to be acceptable in the construction of a model and 
suffices for the needs of Sec. 6. This restriction could likely be relaxed with sufficient 
effort. 

Let » denote the measure generated by the non-decreasing function — In M of 
an admissible J/, 

An admissible trajectory is any ordered pair of real functions x, y of ¢ defined by an 
admissible M and 6 and the system of integral equations 


“0 


nt ) 
u(t) =c | cos Ody + u, | 
| 


at 


v(t) =c | sin 6 du + v, 
7 (2.1) 


at at 
z(t) = | u(r) dr, y(t) = v(r) dr. 

0 “0 4 
If M should have a continuous derivative J/’ then dy could be replaced by —M’ dr/M 
in (2.1). Familiar trajectory equations ((6, p. 172]) would follow by differentiation. 
Constant c is to be interpreted as magnitude of the exit velocity. Functions u(é), 
v(t) introduced is (2.1) are discontinuous when and only when M is discontinuous. Clearly 


u(0O—) = u, v(O—) = v. Moreover | u(t) |, | v(t) | are bounded by c In (W(O—)/M(T+)); 
hence x(t), y(t) are Lipschitzian and therefore AC on any time-interval. We then know 
that x and y are differentiable a.e. (almost everywhere) and that 2x’(t) u(t), y(t) = 


v(t) a.e. on [0, 7] or any other interval. 
If M, @ and the pair z, y are all admissible the ordered quadruple (M, @, xz, y) will 


be termed admissible. 
ProBLEM. To minimize M(O—)/M(T+) on the class of all admissible quadruples 


satisfying the fixed terminal conditions 
u(T) = U, oT) = V, z(T) =a, y(T) = b. (2.2) 


Both the existence and the nature of a minimizing quadruple will be established 
via an auxiliary problem. 

3. Formulation of the auxiliary problem. Let n > 2 be a fixed positive integer 
and consider the new problem obtained from that of Sec. 2 by restricting M to be a 
step-function with at most (nm + 1) steps at 4) = 0,4, = 7’, and at times t, ,t,, --- , ta-1, 
which may be anywhere on the closed interval [0, 7] and are not necessarily in the order 
indicated by subscripts. Let r; denote the mass-ratio at ¢; , viz. 


r, = M(t; —)/M(t,; +), t= 0,-:- ,n. 
Integrals (2.1) are now expressible as sums on 7 from 0 to n and the auxiliary problem 


can be rephrased as the following ordinary minimum problem. 
AUXILIARY PRoBLEM. T'o minimize 


> Inr,; = In (M(O—)/M(T+)) (3.1) 
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subject to the constraint 
bs (T — t,)(cos @;) nr; = (a — uT)/c,] 
> (7 — t)(sin 6) Inr; = (b — vT)/e, 





> (cos 6,) Inr; = (U — u)/e, (3.2) 
> (sin 6,) nr, = (V — »)/e, 
and to the further restrictions 
m2 i, i=0,1,-:-,n, (3.3) 
i = 0, tj=T and 0<:t, <T, t¢=1,---,@-— }). (3.4) 


Clearly the infimum of M(0—)/M(T+) does not excede the particular value K 
obtained by setting r; = 1 and assigning arbitrary values to 6; ,7 = 1, --- , (n — 1), 
and completing the elementary solution of the four equations (3.2) for ro , Ta , 8 , % - 
Neither the existence nor nature of solutions of the auxiliary problem is affected by 
adjoining the restriction 


, & My 7=0,:--- ,n. (3.5) 


Since @; enters only through the periodic functions sine and cosine we are similarly 
free to add the restriction 


0 < 6; < 2z, ¢=Q, --+ an. (3.6) 
Side conditions (3.2) through (3.6) define a bounded closed non-empty subset G 

of (3n + 1) — space; hence there exists a minimizing point 
(To a eli Ts ; A ; 6, ee a = 4 ; t oS ~~ s (3.7) 


in G for the continuous function = In 7; . 

Successive reductions of the minimizing set of values (3.7) are effected as follows. 

If t; = 0 drop the jth term in each equation (3.1), (3.2) at the same time replacing 
In ro in (3.1) by In ror; and adding jth terms to Oth terms vectorially in (3.2). Side con- 
ditions (3.2), (3.3), (3.4) continue to hold but the number of summands in (3.1), (3.2) 
is reduced by unity. Similar remarks apply if t; = 7 or more generally if ¢; = & . 

If r; = 1, 7 + 0 or n the jth terms in (3.1), (3.2) are understood to be deleted. It 
is convenient to retain first and last terms in these sums as objects of discussion even 
when such terms equal zero. 

Finally if @; = 6, mod 2x, j7 + k, replace the pairs of terms in (3.1), (3.2) by a single 
term in each case using the solution ¢, r of the system 


(T—dlnr = (T — ¢t) Int, + (T — 2, Inn, , 


Inr = Inr; +Inr,. 


One verifies that 
(¢—t):(4 -—0 = Inn: Inr,; = ( —):(t- 


hence that ¢ is on the interval [0, 7’). 
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Reductions of the several types can be applied and remaining values (3.7) relabeled 
so as to yield a reduced set 


Vesta 5?" Tm» GO, °° » On yt, °° » Meer msn, (3.8) 


satisfying the restrictions 
fo, = I, fi 3: i> 4, t=1,---,(m—-)), (3.9) 


6; ~ 6 mod 2r, j#k, (3.10) 


| 
St, < Ff, +=1,°---,(m—- l), | (3.11) 
t~At, jk. | 


By the reduced auxiliary problem or reduced problem for short we mean the problem 
of minimizing (3.1) subject to side conditions (3.2), (3.9), (3.10), (3.11). The range of 
summation will now be denoted by m as a reminder of the reduction process. Since 
the auxiliary problem has a solution the reduced problem necessarily has a solution. 
Moreover a solution of the reduced problem is a solution of the original auxiliary problem. 

4. Application of Lagrange multipliers. Consider the problem f(z) = minimum 
subject to constraints ¢,(r) = 0, --- , g.(x) = 0, in which x denotes (x, , --- , 2) 
and a < @. The usual multiplier process yields a necessary condition on a minimizing 
point X provided that (see for example [2, pp. 544-547]) X is an interior point of the 
set on which a minimum is sought and provided that the a by 8 matrix of derivatives 
of ¢; , *** , ¢ With respect to x, , --+ , 23 is of rank a. 

In the proofs of Sec. 5 the interior nature of the minimizing point follows from 
selected strict inequalities (3.9), (3.11). One can verify that the matrix has the required 
rank in each case. Further mention of this condition omitted. 

Let R, A, B, C, D denote respective left members of (3.1), (3.2) for the reduced 
problem, let X, , - , A, be multipliers, and let F = R+A,A + A.B + A.C + AD. 


We record for reference that 


dF /dr; = — [1 + A,(T — ¢,) cos 6; + A(T — ¢,) sin 6; + Az cos 0; + Ay sin 4,], 
7 


dF /d0; = Inr,[—A,(T — ¢,) sin 6; + A.(T — ¢,) cos 6; — As sin 0; + A, cos 8]. 


dF /dt; = —Inr,[dA, cos 6; + Az sin 6;]. 


We shall use certain of the equations 


oF /dr; = 0, p= G, «+=. om: (4.1) 
OF /d0; = 0, gat 0, «+s 2m; (4.2) 
oF /dt; = 0, t=1,--- ,(m— 1). (4.3) 


5. Properties of a solution of the reduced auxiliary problem. 


Let 
Se = ae a Ses oe (5.1) 


denote a solution of the reduced problem. 
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THEOREM 5.1. A solution (5.1) can include at most one mass ratio r; greater than unity 
whose subscript differs from 0 or m. 

Proof. Suppose that r; > 1,7, > 1,7 # k, and that neither j nor k is 0 or m. Con- 
sider the minimum of (3.1) in six variables r; , 6; , t; , 7 = j, k subject to four constraints 
(3.2) with other variables fixed at values (5.1). 

From (4.3) with i = j,k 


0, 
0. 


A, cos 0; + A, sin 8; 


Ay cos 6, a dz sin 6, 
It follows that either 4, = A. = 0 or in view of (3.10) that 6; , 6, differ by an odd 


multiple of 7. 
Under the first alternative we see from (4.1), (4.2) with 7 = j that A; = —cos 0; 


and \, = —sin 6; . Similarly with 7 = k, \; = —cos 6, and \, = —sin 6, . Such results 
contradict (3.10); hence cannot occur. 

Under the second alternative we have from (4.3) and (4.1) for 7 = j that A, cos 
6; + A, sin 8 —1 and for 7 = k the contradiction that A; cos 6, + A, sin 6, = —Ag; cos 
6; — 4 sin 8; —1. 

TureoreM 5.2. Jf ry > 1 ina solution (5.1) thenr; = 1,4 = 1, -++ , (m — 1). 


Proof. Suppose that r, > 1 and that r; > 1 for some fixed j between 0 and (m — 1) 
inclusive. Consider the minimum problem in five variables ro , r; , % , 9; , t; with all 
other variables fixed at values (5.1). From the five equations consisting of (4.1), (4.2) 


for 7 = 0, j and (4.3) for 7 = j we obtain relations 

A, cos 6; + A, sin 8; = 0, (5.2) 
A; cos 6; + Aysin 6; = —1, (5.3) 
T — t)(A, Cos 0 + A. Sin 8) + A; Cos O% + A,sIn & = —1, (5.4) 
(T — t,)(—X, sin 0 + Az Cos O) — Az SiN OH + Ay Cos 6 = O, (5.5) 
(T — t,(—X, sin 0; + A. cos 6;) — Az sin 6; + A, cos 6; = O. (5.6) 

Solving (5.2) and (5.5) for A, , A. we find that 
d —(sin 6,)(Az3 Sin 8) — Ax COs O)/(7’ — t,) cos (6; — 60) (5.7 
4. = (cos 8,)(Az3 SIN A — Ay COS O)/(T — to) cos (0; — 6,).J . 
If cos (8; — 0) were zero then d; sin 0 — d, cos 6 = 0, which contradicts the result 


obtained from (5.3) by replacing cos 6; , sin 6; with —sin 4) , cos 6 . Thus denominators 
in (5.7) cannot vanish. 
From (5.6), (5.7) 


Az [sin 6; — (7 — t,)(sin 6)/(T — to) cos (0; — 4%)] (5.8) 
+ r,[—cos 0; + (T — t,;)(cos 0)/(T — to) cos (0; — %)] = 0, 
while from (5.8), (5.3) with the fact that t = 0, 
A; = (T/t;))[—cos 0; + (T — t,)(cos 6)/T cos (0; — %)], (5.9) 
A, = (7/t)[—sin 6; + (T — t,)(sin 6)/T cos (6; — 4)]. 
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Substitute (5.9) into (5.7) and find that 
\4, = —?¢;'sin 6; tan (0; — 6), \. = t;' cos 6; tan (0; — 6). (5.10) 


Finally substitute (5.10) and (5.9) into (5.4), which reduces to cos (0; — 6) = 1, 
in contradiction with (3.10). 

THEOREM 5.3. Jfr, > 1 ina solution (5.1) thenr; = 1,4 = 1, --- , (m — 1). 

Proof. The proof follows that of Theorem 5.2 as far as (5.6) with subscript 0 re- 
placed by m. We now understand this change to have been made in (5.2), --- , (5.6). 

The proof now diverges as a result of the fact that (7’ — t,,) = 0. From (5.4), (5.5) 


accordingly reduced \, = —cos 6, , A, = —sin 6, , which with (5.3) again contradicts 
(3.10). 

The following theorem summarizes results of Sec. 3, 5. For completeness the trivial 
case in which 7 = 0 and m = 1, thus far excluded from our formulation, is included 
as (v). 

THEOREM 5.4. The auxiliary problem necessarily has a solution of one of the following 
types 


(i) ro > 1,7, > 1, each other ratio = 1, 

(ii) 1 > 1, each other ratio = 1, 

(iii) r,, > 1, each other ratio = 1, 

(iv) r; > 1 for one 7 ¥ 0, m and each other ratio = 1, 

(v) m = 1, ro > 1 and each other ratio = 1. 

THEOREM 5.5. Jf the auxiliary problem has a type (iv) minimum it also has a type 


(i) minimum. 
Proof. Let h, H respectively denote the expressions av — bu, aV — bU and let 


R denote the positive square root of (U — u)* + (V — v)’. If a type (iv) solution exists 
elementary relations among distance, speed, and time imply that 
T = (h — H)/WU — uV). (5.11) 
Moreover 


c (min bd In¢y;) = #. (5.12) 


If we now solve (3.2) with m = 2 and r, = 1, taking the preassigned 7’ to be (5.11) 


we find that 

ec > Inr, = R{[W’/(h — HY]? + [H’/(h — H)’]'”}. 
The stated conclusion follows from (5.12) and the fact that the above expression in 
braces cannot exceed unity. 

6. Solution of the original problem. Let (M, 0, x, y) be any admissible quadruple 
satisfying terminal conditions (2.2). It is convenient to the application below of a 
convergence theorem for Stieltjes integrals to consider a sequence M, of step functions 
converging to M or the interval [—1, 7 + 1] except possibly at discontinuities of M. 
(Any other closed interval to which [0, 7] is interior would do.) Such a sequence M, 
can be obtained using a suitable nested sequence P, of partitions of [0, 7]. Define M,(t) 
for —1 < ¢t < Oas M(0O-—), for T < t < T + 1 as M(T+) and for each of the disjoint 
half-open subintervals whose union is the half-open interval [0, 7] as the mean of M 
over that subinterval. Let 6, be any sequence of continuous functions converging uni- 
formly on [—1, 7 + 1] to 6; in particular take 0, = 0,v = 1,2, ---. 

A well-known theorem ((3, p. 283, Theorem 23]) insures that the sequences 
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t 
u,(t) = c¢ [ cos 6, du, + u, 
“0 (6.1) 
t 
c [ sin 6, du, + v, 


“0 


v,(t) 


converge for te [—1, 7’ + 1] except for / a point of discontinuity of M. (See Hypothesis 
B, , [3], p. 281). Moreover limits of right members of (6.1) are right members of the 
first two equations (2.1) written for the M, @ given above except for the countable 
set of discontinuities of M. 

The ordinary bounded convergence theorem then insures that 


2if = i u(ndr, ()= [ v,(1) dr, (6.2) 


converge to the admissible trajectory x, y included in the admissible quadruple with 
which we started. 

The position (a, , b,) given by (6.2) and velocity (U, , V,) given by (6.1) fort = T 
may differ from values (2.2) but u,(7'+), v,(T7+), z,(T), y,(T) necessarily converge 
to respective values (2.2). 

By Theorem 5.4 there is a step function //* with at most two steps and a continuous 
function 6* whose values at the discontinuities of M* together with M* define a solution 
to the auxiliary problem with terminal position (a, , b,) and terminal velocity (U, , V,) 
at 7. Hence 

In [M*(0—)/M*(T+)] < In [M,0—)/M,(T+)]. (6.3) 

Each /* is of one of the types enumerated in Theorem 5.4; hence through possible 
extraction of subsequences we can suppose all sequences on v so chosen that each M% , 
vy = 1, 2, --- , is of the same type and such that M* converges to a limit M* on [—1, 
T + 1], which is necessarily a step function with at most the same number of steps 
as M* . Function 6* can be taken as constant or linear according as M* has one or 
two steps. We can suppose 6* convergent to a limit 6*, again by restriction to suitable 
subsequences. Terminal position x*(7), y*(7) and velocity u*(7+), v*(T7+) deter- 
mined by M* , 6* , and (2.1) converge to values (2.2) as observed above. 

Passage to the limit in (6.3) shows that 


In [M*(O—)/M*(T+)] < In [M(O—)/M(T+)]; 
hence that the particular admissible quadruple (M*, 6*, x*, y*) is at least as good as 


(M, 6, x, y). 

Accordingly the minimum problem stated in Sec. 2 has as a solution any solution 
of the auxiliary minimum problem. 

To solve an example refer to the possibilities of Theorem 5.4. Since at most four 
quantities are subject to four side conditions (3.2) multipliers are not applicable at 


this stage. One or more solutions to the example are provided by (3.2) with n = 2. 
7. Examples. 


(u, v) (U, V) (a, b) T type c (Min = In r,) 
(0, 1) (0, —2) (1, 0) l (i) 7” +5” 
(1,0) (2,0) (2, 0) 1 (ii) 1 

(1, 0) (2, 0) (2, 0) 2 (iii) 1 

(1, 1) (2, —2) (2, 0) 2/3 (iv) 10°”? 

(1, 1) (4, —2) (2, 0) 1 (i) or (iv) 3 (2)'” 
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8. Variable transit time. The minimum of M(0—)/M(7T+) denoted by dX is a 
function of transit time 7’. If the problem is altered by allowing 7' to vary over the 
non-negaitve reals the solution when it exists is available by the present paper plus 
@ minimization with respect to 7. There are cases in which no minimum exists, e.g. 
(u,v) = (0,0), (U, V) = (-—1, 0), (a, 6) = (1, 0). For each fixed 7 we have a Type 1 
solution. Inf \(7’) is approached as 7 — © but is never attained. 

9. Concluding remarks. The reviewer has called attention to a forthcoming paper 
by E. W. Graham in the Journal of the Aero/Space Sciences which may overlap the 
present results and which includes uniform fields. That the uniform field case is mathe- 
matically equivalent to the free space problem can be seen as follows. In the first two 
equations (2.1) add respective constant multiples 2A7, 2BT of t; hence in the last 
equations (2.1) add AJ’, BT’. In the respective equations (3.2) we must then replace 
a, b by a(l — T”’), b(1 — T”) and replace U, V by U — 2AT and V — 2BT. Subject 
only to these substitutions all present results apply. 

Among questions for which only partial theories appear to exist are the following. 

If terms representing a general gravitational field are included in (2.1) as in [4], 
for what non-uniform fields is there an analog of Theorem 5.4, eliminating the need 
for maneuvers with more than two null-thrust arcs like that of [6, p. 175, Fig. 6]? 

If our problem is modified into Type 2 by adjoining a Lipschitz condition on M, 
is an optimal maneuver realizable with at most two periods of maximal thrust? 

With aerodynamic forces also present we have the possibility of time intervals on 
which M decreases continuously in addition to possible coasting periods together with 
impulsive burnings or periods of maximal thrust according as the problem is of Type 1 
or 2. An extension of Theorem 5.4 coupled with proof of the existence of a best program 
should be of interest. 
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NON-DECREASING STEP RESPONSES WHOSE TRANSFER 
FUNCTIONS ARE SUBCLASS K * 


BY 
ARMEN H. ZEMANIAN 


New York University 


Introduction. ‘This paper is concerned with the intersection of two classes of rational 
system functions. The first class was called the subclass k in previous papers and its 
functions constitute generalizations of the positive-real functions [1-4]. Functions of 
this type are quite common in electric network theory. For instance, the reciprocals of 
Hurwitz polynomials are members of this class and they are used quite frequently in 
network synthesis [5-7]. The transfer functions of the series-shunt peaking network 
and the Dietzold network [8] are also subclass k functions. Or again, a minimum-phase 
function having a monotonic decreasing** phase function is a subclass k function [1]. 
The unit step responses of such functions need not be monotonic. 

The second class is comprised of those rational system functions whose corresponding 
unit step responses are non-decreasing [9-12]. Such functions are of interest in those 
applications where an overshoot in the step response is to be avoided. For instance, this 
is often the case in measuring apparatus. 

It has been shown previously that certain restrictions exist on the step response of 
any subclass / function. If, in addition, this response is known to be monotonic increasing, 
it is possible to strengthen these restrictions. This is one of the objectives of this paper. 
Of course, the results that will be developed here are useful only if it is known that the 
step response is monotonic increasing. Various sets of sufficient conditions that insure 
this are given in [9-12]. 

The second objective is the following. In [1], lower bounds on the rise time and the 
settling time of the unit step response corresponding to a subclass k transfer function 
were obtained. These lower bounds were the envelopes of families of straight lines 
wherein the coefficients of the straight line expressions contained certain infinite series. 
Through numerical computation these series were summed and then the envelopes were 
determined by graphical means. Results were obtained for values of k that ranged from 
one to five. In this paper the infinite series are summed into closed forms and parametric 
equations for the envelopes, which hold for all positive integer values of k, are developed. 
Such expressions for non-decreasing step responses are given within the body of the 
paper. The parametric equations for the envelopes when the step responses need not be 
monotonic (i.e., such as those of [1]) are given in the appendix. 

Finally, having such results available, one may formulate lower bounds (in terms of 
definite integrals of certain functions) on the moments of non-decreasing step responses 
having subclass & transfer functions. This additional result is also of significance in 
probability theory wherein moments have been used to describe the probability dis- 
tribution function. In this case a distribution function corresponds to the non-decreasing 


* Received February 8, 1960. 
**Here, the phase angle is considered positive when a steady-state response leads a sinusoidal input. 
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Fig. 1. Illustration of the settling time 73 (a) for a unit step response that has overshoot and (b) for a 
unit step response that is non-decreasing. 


step response and its characteristic function corresponds to the system function. The 
rationality restriction on the characteristic function or on the system function may be 
dropped if the subclass k functions are defined as in [3]. 

The notations used here will be the same as those of [1] and this paper should be 
read as a sequel to it. 

Definitions. Let Z(s) be the following rational system function of the complex 
variable s = o + jw, where all the coefficients are real numbers 


s+ a8"? +--+ +a _ ~N@ (1) 


s” + bys” + ees $b Ds) 

Furthermore, let K > 0,k = m — n> 1, Z(0) = Ka,/b, = r, and let D(s) be a Hurwitz 
polynomial. The real and imaginary parts of Z(jw) will be denoted by R(w) and /(w), 
respectively. Under these conditions, the response A(t) of the initially quiescent system 
to a unit step input applied at time t = 0 will be a continuous function for all ¢. A(t) is 
related to Z(s) through the following Laplace-Stieltjes transform 





Z(s) = 


a) 


As) = | ddd. (2) 
70 

The subclass k is defined as follows. Z(s) will be called a subclass k function if the 
following conditions are satisfied. For k odd, R(w) has k — 1 changes of signin —- < 
w < © and R(w) is positive in the neighborhood of w = 0. For k even, I(w) has k — 1 

changes of signin —©- <w< @ and dI/dw is negative in the neighborhood of w = 0. 
In the following analysis we shall assume that the final value r of the unit step 
response is a positive quantity. Another quantity that will be considered is the settling 
time to 6 which is denoted by 7; . This settling time is the least time beyond which A (¢) 
remains within the bounds (1 + 6)r where 0 < 6 < 1. A particular settling time is 
illustrated in Fig. 1 for a unit step response that has overshoot and for another unit 
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1.0 





Fig. 2. The lower bounds on the settling times for the unit step responses of the subclass k functions. 
The solid curves are the lower bounds in the case where the unit step responses are non-decreasing. 
The dotted curves are for the general case where the unit step responses need not be monotonic. 


step response that is non-decreasing. When A(t) is non-decreasing, any value of time 
t is a settling time 7; where 6 = 1 — A(t)/r. 

Upper bounds on non-decreasing step responses. If a system function is a subclass 
k: function, the settling time to 6 cannot be less than certain lower bounds given by the 
dotted curves of Fig. 2 [1]. That is, given the unit step response of a subclass k function, 
any deviation dr (0 < 6 < 1) from the final value r will determine a settling time 7; ; the 
point, whose coordinates are 7;/7'y and (1 — 4), must lie to the right of the dotted curve 
for the given k. In Fig. 2 the time scale is normalized by the factor 7'y which is given by 


n= [x] @ 


For non-decreasing step responses this result can be strengthened. The new lower 
bounds on the settling time are given by the solid curves of Fig. 2. They may be derived 
in a fashion similar to the development of the dotted curves. In fact, the first parts of 
the two derivations are precisely the same and will not be repeated in this paper. This 
section gives the derivation of the solid curves starting at the point where the develop- 
ment diverges from that of the dotted curves. 

This result and the monotonicity of A(¢) together imply a still stronger restriction. 
Any time ¢ is a settling time 7; for non-decreasing step responses. Hence, for a given k 
the corresponding solid curve of Fig. 2 is an upper bound on a non-decreasing A(t). The 
upper bounds that were previously shown to exist [1] on the entire unit step responses 
of all subclass & functions are given by the dotted curves of Fig. 3. For comparison, the 
new upper bounds are repeated in this figure as the solid curves and a substantial improve- 
ment may be noted. The examples of the next section indicate that these upper bounds 
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Fic. 3. The upper bounds on the unit step responses of the subclass k functions. The solid curves are 
the upper bounds in the case where the unit step responses are non-decreasing. The dotted curves are 
for the general case where the unit step responses need not be monotonic. 


are fairly close to being best possible when k = 1 except possibly in the neighborhood 


of the final value line. 
This conclusion may be stated as follows. 
Theorem. Let Z(s) be a subclass k function, let the corresponding A(t) be non-decreasing, 


let B be the upper bound on A(t)/r, and let 


| t |* 
Fea —|——i. (4) 
* ik! | 


Then the parametric equations for the upper bound on A(t) are given by (4), (5), and (6) 
for k odd and by (4), (7), and (8) for k even, where 0 < y <1. 


(k-3)/2 —. 1 k i 
t= 2 "1-7-3... } 1 - -—— uy e-2p( 1) 
2 > y | - al 1 — ry ctn ry tan") ir) 
(> 


a 

















y as bi 
+ 30 — ay cin xy) (* TULA) — 4 log 2] + 2')}, 
Ba —— te tig SF" eG ayy ts ~ be 
~ rl — wy ctn wy) - —_— = “H)Y ok 2u-1 Ve-2y . 
)) 
QA yO’ (y) 
— 2log2+ i) 4 wr) 
er i k-2 ] k —_ Qu ° \ - 
3 = — 2 ( es saiee: 2 er = Ss spin —? ( 9 { 
f=14 dX y E po aE ae = i »,(1) 7) 
2 sin TY Om”? k-2 “| ] | ° 
B=1- — (k= * oo mn Bef lh (8) 
m1 — ry ctn ry) 2d e — Su)y i ot-2e-1 te-26(1), 
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where 


2%) = v(1 Z u) v(2 - u) + oi 4 u) ~ v(4 + u) : (9) 


vw =3[-w(0-H)+¥G-s) +e +8) -9G+))] oo 


-_ 


In (5) and (6) the summations on u do not exist if k = 1. Furthermore, Li,(1) is the poly- 
logarithm of order p and argument one and V(x) and ¥'(x) are the digamma function and 
its derivative, respectively. 

Note. Tables of Li,(x) and of V(x) and W(x) may be found in [13] and [14], respectively. 
It should also be noted that this theorem holds for a wider class of transfer functions 
i.e., for those functions that satisfy inequalities (8) and (9) in [1}). 

Proof. As mentioned previously, the first part of the proof for k > 2 is exactly as 
given in [1] up to expression (35) and this part will not be repeated here. For k = 1, see 
[15]. Expression (35) of [1] is the following inequality wherein 0 < y < 1. (Note that 
the upper limit NV in the summation of expression (35) of [1] should be replaced by ~.) 


A(t) < sin ry [xe +} Qy**? D2 (<1 Aya) |. (11) 
my Lk! - fe Pe-v) 

Now those values for A(pt/y) that maximize the right-hand side of (11) must be found. 

Since 0 < y < 1, (sin ry)/ry is positive and the coefficients of the A(pt/y) alternate in 

sign and decrease in magnitude, the first coefficient being positive. Furthermore, the 

A(pt/y) are non-decreasing and no larger than r. Hence, for any given, permissible set 

of values for A(vt/y), (v = 1, 3, 5, ---), the summation will assume its maximum value 


when A[(v + 1)t/y] is set equal to A(vt/y). Under this condition the summation in (11) 


may be rewritten as 


zt) ] 1 } 
2y ° A ny — ae a. ee eae ke ee 12 
> (! — —y) (y+ 1)'[(~@+ 1° -— y] (12) 
» odd 

The coefficients of A(vt/y) in (12) are all positive so that this expression and hence the 
right-hand side of (11) are maximized by setting A(vt/y) = r. 

Inserting the quantities B and £, which are defined in the hypothesis, and rearrang- 
ing the result, (12) becomes the following equality. 








i _ Ty _. Gaih*9 — _(-17"""_ 
t= aie ry B “y 2 pp aah y’) (13) 
Let 
asi ry 
( = — 
(y) sin ry (14) 
and let 
’ +2 = (—1)?*? 
E(a = 2 P “kk; 2 ay 15 
ait 27 -y) (15) 


so that (13) becomes 


& = C(y)B — Ey). (16) 
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Now E(y) may be summed by first expanding each term of the series into partial fractions. 


1 13 FE = | 1 | 1 - £0"). 
+ = 1 = i (17 
pip—y) 4% U+ (“0 ‘yt + oy lp—y pty 7) 


Summing the last two partial fractions in (17) according to (15), we have from [16; p. 
20, formula 1.7.4.(33)] 











— p+1 me a _ a —_ 2(y) 4 
2 (=1) [2 +-1.|- 2 oe 


p=1 
= : l 1 
SE) tid mn 
2, ( : [—- | 
=. y 1 y y re 
=A fo( - ¥) 93-4) - v1 +H) +u(b+4)] 


The right-hand side of the last equation may be simplified by using [16; p. 16, formulas 


1.7.1.(11)] 
tan! + etn | — 1 . aoe § (19) 
2 2 y sin ry y 


a: ee 
2d (- ye |- 


p=1 








dla 


Moreover, 


= (—1)”*' : 
>» oar — = log2 (20) 
p=1 ) 

and, for k = 2, 3,4, ---, 


— 
We dd-2 doe 


p=1 p=1 DP =1 (2 


The right-hand side of this sina may be summed according to [13; p. 169, Eq. (7.1)] 
(— ] ce 


: —_ = [ _ sh |uis). (21) 
p=) Dp ys 


Combining (15) and (17) through (21), a closed form for E(y) is found. For k odd 
2¢(y) 


E(y) = -—2 2. y* “ly = ai [ea ou(1) — 2y log 2 + ys > (22) 


and for k even 


(k-2)/2 
4 ( | = —92 a? —_- Li. ~~ + —— — 23) 
E(y) > y E ah be (1) 4 tues f. (23) 





In (22) the summation on yu does not exist if k = 1. 
The parametric equations for the desired envelopes may be obtained by using a 
standard procedure [17; p. 85]. Differentiating (16) with respect to y, one obtains 
E'(y) 
— (24) 
Cy) 


Differentiating (14), (22), and (23) and combining them according to (24), we obtain 
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(6) and (8). Furthermore, inserting (24) into (16), we have 


_ C(yk'Yy) 


Through some elementary manipulations, (25) may be converted into (5) and (7). 
This completes the proof. 

Before leaving this section, it should be noted that in certain cases the upper bounds 
on a non-monotonic A(t) given by the dotted curves of Fig. 3 can also be strengthened. 
This is possible if it is known that A(é) satisfies the following conditions. 


(a) A(t) is non-decreasing in the interval, 0 < t < 7, , where A(r,) <r. 
(b) Beyond 7, , A(t) remains within the bounds, 


A(r,) < A(t) < 2r — A(r,). 


In this case the same reasoning used before leads to the following conclusion. In the 
interval, 0 < t < +, , A(t) is bounded according to that dotted curve of Fig. 2 which 
has the appropriate k. Unfortunately, for a given Z(s) a value of 7, that satisfies con- 
ditions (a) and (b) is not usually known so that this last result is not readily applicable. 

Examples. ‘To get an estimate of how close the bounds of the theorem are to being 
best possible, two examples, were computed for the unit step responses of the subclass 
| functions. The first subclass 1 function is the normalized driving-point impedance of 
the critically damped, shunt-peaked filter, given by the following expression 


s+4 
Z(s) = oo (26) 
The corresponding unit step response is 
A(t) =1—(1+ de. (27) 
It is plotted as curve (b) in Fig. 4. Curve (a) of Fig. 4 is the upper bound on such unit 
step responses; i.e., it is the solid curve for k = 1 in Fig. 3. Any possible improvement 


in this bound cannot be greater than the amount indicated by curve (b). 

It can be noted that the horizontal distance between curves (a) and (b) is greatest 
in the neighborhood of the final value line. Another example of a non-decreasing unit 
step response, which is closer to the upper bound in this neighborhood and corresponds 
to a subelass 1 function, is provided by the following 


(s + 2.56)” + (3.56)° 











‘2 oS ere te 28 
4(s) (s + 2.56)[(s + 2.56)" + 1] (28) 
The corresponding unit step response is non-decreasing and is given by 

A(t) = 1 —e* [4.98 + 1.56(sin ¢ — 2.56 cos 2]. (29) 


It is plotted as curve (c) in Fig. 4 and provides a sharper estimate on any possible improve- 
ment in curve (a) for A(t)/r > 0.9. 

The moments of non-decreasing unit step responses of subclass k functions. The 
moments of non-decreasing unit step responses are given by 


i. 
m=—[ tA, 
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Fia. 4, (a) The upper bound on the monotonic unit step responses of the subclass 1 functions. (b) The 
unit step response of the shunt-peaked filter given by Eq. (27). (c) The unit step response given by 


Eq. (29). 


0 < t < ~, lower bounds on m, may be obtained by replacing A(t)/r by its non-decreas- 
ing upper bound B(y). Note that B(1) = A()/r = 1. Furthermore, using (4) we have 


where r = Z(0) = A(~) andi = 1, 2, 3, --- . Since ¢’ and A(t) are non-decreasing in 


al 


oil 
m. > Ti | (ke) 4 om 


v0 
Z(s) is rational so that A(t)/r cannot achieve its upper bound. Hence, equality cannot 
be achieved in this expression. Here, &(y) is given by either (5) or (7) and dB/dy is the 
derivative of either (6) or (8). Thus, integral expressions for the lower bounds on these 
movements have been formulated. 

APPENDIX 


PARAMETRIC EQUATIONS FOR ENVELOPES DEVELOPED IN [1] FOR THE CASE 
Wuere A(t) Is Not NEcEessartty MONOTONIC 


1. The lower bounds on the rise time. As in [1], let us define the rise time 7’, as 
the least time where A(7',) = r > O and beyond which A(t) is bounded by (1 + y)r, 
y being a positive number. Expression (36) of [1] gives a family of straight lines which 
define an envelope that yields a lower bound on 7’, . This expression is repeated here. 
For 0 < y < 1, we have 


.if.f : 
kl [| > —Qy)v + Vy), (30) 
where 
vs k+2 — — 
Qty) = 2y Lig — yy’ (31) 
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7 — —- k+2 = = (—1)"** 
b (y) = ain ry 2 pr rip fon y’) (32) 


Let £, be the lower bound for the left-hand side of (30). Thus, 


7 


4 a (k!¢,)'”* (33) 
N 
and (30) becomes 


E, = —Qiy)y + V(y). (34) 


Proceeding as in [17; p. 85], the parametric equations for the envelope of the family 
of straight lines of (34) are found to be 





_ vy) 
"Tw ™ 
—- VY , y . 
t, = Ga? + Va. (36) 


The expressions for Q(y) and its derivative may be summed into the following closed 
forms by using a procedure similar to that employed in summing F(y). For k odd, 


(k-3)/2 
Qy)=-2 Do y* *Lix-2,(1) — y[W(y) + W(—y) + 2g], (37) 
k-3)/2 
Q’(_y) = -2 D> (k — Quy Liz-2,(1) — Vy) 


(38) 
— W(—y) — 2g — yW'(y) + y¥'(-y), 


where g is Euler’s constant (0.5772 ---) and where the summation on u does not appear 
if | 1. For k even, 





(k-2)/2 
Qy) = —-2 do y* *Liz-2,(1) + 1 — ry ctn zy, (39) 
u=0 
QO" 4) = —2 5 Bho 2 ; ieee a ° _ . = = Ss | d 
I(4 2 a (h uy Lix-2,(1) + no E = cos ry (40) 
Furthermore, by comparing (32) with (14) and (15) we see that 
V(y) = Cly) — Ey). (41) 
Hence tor } odd 
Viy=2 dX yr*1-- = Li,-2,(1) + 2y log 2 — yy) + =. = (42) 
pre 2, | gk-2u-1 | b— 2p “y = 3 sin ry ’ 
> l ; 
V'(y) = 2 a k— Quy" [1 ae [hiss 
; , (43) 
Q( Q 
+ 2 log 2 — ay) _ y2'y) + —"— (1 — zy ctn ry), 


2 2 sin ry 


where 2(y) and Q’(y) are given by (9) and (10), respectively, and where the summation 
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on uw does not exist if k = 1. For k even 
k—-2)/2 2 1 
V /) as l + 2 > y an = 9* Lois, (44) 


- > ae ee l = 
V'"(y) = 2 z. (k — 2Qy)y [1 ~ = Joana. (45) 


u=0 


Thus, (33) and (35) are the parametric equations for the lower bounds on the rise time 
where £, is given by (36). The quantities in (35) and (36) are given by (37) through (45), 
where 0 < y < 1. For & ranging from one to five, these bounds are plotted in Fig. 3 
of [1]. 

2. The lower bounds on the settling time. ‘The parametric equations for the lower 
bounds on the settling time may be developed in a similar fashion. From expression 


(39) of [1], we have for 0 < y < 1 


l T§ ' , 
: = : i Pi 
kl | > —J(y) 6+ Vy), (46 
where V(y) is given by (32) and 
TY M2 € l ¥ 
J(y) = —* +2" > ; = (47) 
Sin ry =P — yy) 


Letting &; be the lower bound for the left-hand side of (46), the parametric equations 


Ss 


for the lower bounds on the settling time are found to be 


p> (ke iar (48) 
N 
V’(y) 
= y a (49) 
e Ly 
where 
V'(y) J (y) ie 
t= —— Hs Wy). (5 
- J ‘(y) r Vy) a 


The closed form expressions for V(y) and V’(y) are given by (42) through (45). Further- 


more, 
J(y) = me + Qy), (51) 
sin wry 
; T > - 
Jy) = —— (1 — ry ctn ry) + Q’(y), (52) 
sin ry 


so that the closed form expression for J(y) and J’(y) are obtained by combining (51) 
and (52) with (37) through (40). For & ranging from one to five, these bounds are plotted 
in Fig. 4 of [1]. 

Acknowledgment. The author is indebted to Professors A. Erdélyi and L. Lewin 
for pointing out how the infinite series in (15) and (31) may be summed. 
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BOOK REVIEWS 


(Continued from p. 354) 


of some classical problems when the stress-strain relations are slightly nonlinear, and nonlinearities 
arising from boundary conditions. 

The subject matter is interspersed with interesting and novel examples, and all ‘‘analytical,’’ numeri- 
cal and graphical methods which are described are illustrated by examples most of which are original. 

A book review of conventional length can not do justice to so comprehensive a book as that under 
review. Therefore, some general impressions will have to suffice. 

The level of the book is such that a reader familiar with elementary calculus and with ordinary 
differential equations is prepared to read it. In fact, the author has purposely addressed himself to the 
practitioner. Advanced mathematical techniques were avoided by ignoring questions of existence of 
the solutions which are desired, and by assuming that the methods used in their construction will con- 
verge. On the one hand, this approach ignores many results now available; on the other, it is quite 
realistic because the present state of theoretical knowledge of nonlinear differential equations is totally 
inadequate to establish rigorously the validity of so-called analytical solutions that have served so 
well to predict experimentally observed phenomena. 

Many sections of this book contain original work by the author, most of it not published heretofore. 
Moreover, the author has shunned no efforts in presenting a great many very detailed solutions to 
specific problems. While these may make for tedious reading for those more interested in theoretical 
work, they will be welcome additions for the practitioner. 

The author’s presentation follows the German tradition and is not easy to read even for the native 
German; it may be unduly difficult for the American reader. The typography and the illustrations are 
exceptionally appealing and the text seems virtually free from typographical errors which so frequently 
mar first printings. 

The reviewer’s preference would have been to delete entirely the section on nonlinear elasticity. 
For one thing, this section does not deserve the promising title which it bears since it deals only with 
nonlinear stress-strain relations, and for another, the connection between this part and the remainder 
of the book is tenuous. The space gained in this manner could then have been used to present additional 
material now omitted. 

Most notable among these omissions is the use of phase-plane methods for the discussion of the 
transient solutions of non-autonomous systems of one degree of freedom. In addition, the method of 
Poincaré for the discussion of solutions as functions of a parameter might well have been included. 
In fact, some of the examples could have been treated more easily and elegantly by this method than 
was done in the text. 

These remarks are not meant as criticisms. In a lesser book they would have been out of place. 
However, in this book the coverage is so complete, the treatment is so thorough, and considerations 
of space seem to have been so liberal, that the reviewer wishes that he could find all of the important 
material under a single cover. In any case, he can recommend this book very highly as an excellent 
reference work and a welcome addition to the growing literature in the field of nonlinear vibrations. 


R. M. RosENBERG 


Refined iterative methods for computation of the solution and the eigenvalues of self-adjoint 
boundary value problems. By M. Engeli, Th. Ginsburg, H. Rutishauser, E. Stiefel. 
Birkhauser Verlag, Basel, Stuttgart, 1959. 107 pp. $3.95. 

Numerical analysis resembles a physical science in that its theories lose or gain luster in confront- 
ing observations. But in both cases the observational data, which for numerical analysis is concerned 
with accuracy, computational speed, generality of application and, frequently, ease of coding, must 
be skillfully marshaled to be convincing. The first three of these criteria have been applied to a variety 
of computational methods, using a fourth-order plate problem as the principal test, and the authors 
have provided an extremely useful evaluation of the methods examined. 


(Continued on p. 386) 
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SOME VARIATIONAL PRINCIPLES FOR PROBLEMS IN 
HYDRODYNAMIC AND HYDROMAGNETIC STABILITY’ 


BY 
R. C. DiPrma 


Department of Mathematics, Rensselaer Polytechnic Institute 


1. Introduction. A number of problems in hydrodynamic and hydromagnetic 
stability give rise to characteristic value problems with differential equations of high 
order but with constant coefficients. Examples are the sixth order problem that arises 
in the stability of a viscous flow between concentric cylinders rotating in the same 
direction when the gap between the cylinders is small compared to the mean radius 
[1]; the eighth order problem that arises for the same geometry but with an electrically 
conducting fluid and an axial magnetic field [2]; and the sixth order problem that arises 
in the inhibition of convection by a magnetic field [3, 4]. (See [1] and the references cited 
there for a more complete list of such problems.) 

Since the differential equations have constant coefficients the characteristic value 
problems can, in theory, be solved exactly by standard methods. However in practice 
such methods are often laborious and unwieldy. In several cases the characteristic 
value problem can be formulated as a rather unusual variational problem depending 
upon two functions which are related by a lower order differential equation. A variational 
principle of this type was apparently first suggested by Pellew and Southwell [5] in their 
treatment of the classical Bérnard problem—the instability of a layer of fluid heated 
from below. More recently Chandrasekhar has used this method to treat several stability 
problems including those mentioned earlier [1, 2, 3, 4]. Although very satisfactory 
results have been obtained from variational principles of this type the computations can 
become somewhat tedious since it is necessary to integrate a differential equation as 
part of the solution. Further, the variational principles are used in such a manner that 
an exact solution of the characteristic value problem requires the evaluation of a de- 
terminate of infinite order. Thus only approximate (though certainly satisfactory) 
answers are available in practice. : 

In this paper it is shown how these characteristic value problems can be formulated 
as variational problems which can be solved exactly without the necessity of integrating 
any differential equation. The method of solution depends primarily upon the expansion 
of the unknown functions in the variational principles in complete sets of appropriate 
orthogonal functions. It is also shown how these methods can be extended to a non-self 
adjoint boundary value problem by using a variational principle dependent upon both 
the original problem and the adjoint boundary value problem. And in Sec. 5 a method 
of solving the characteristic value problem by a direct series substitution in the differ- 
ential equation is illustrated. 

2. The Taylor problem. ‘The problem of the stability of a viscous fluid between 
two infinitely long concentric rotating cylinders was first successfully investigated both 


1 Received Feb. 25, 1960. This work was sponsored by the Office of Naval Research under Contracts 
Nonr-591(08) and Nonr-591(12). 








376 R. C. DiPrima {[Vol. XVIII, No. 4 


theoretically and experimentally by G. I. Taylor [6]. In the case that the gap between 
the cylinders is small compared to a mean radius, and the cylinders are rotating in the 
same direction, the stability problem reduces to the determination of the minimum real 
positive value of 7 for all real positive a for which the following boundary value problem 
has a solution: 
(D? — a’)*v = —a’To, —3 <2 < 3, (1) 
v= (D?-—a@yv = D(\D—awv=0 at «= +}. (2) 
Here D denotes d/dx, a is a dimensionless wave number, and 7’ is a dimensionless number 
(called the Taylor number) which depends upon the geometry and the angular velocities 
of the cylinders. (Set a = 0 in Eq. 10 of [7] to obtain Eq. (1). See also [1]). The boundary 
ralue problem defined by Eqs. (1) and (2) has been solved by Pellew and Southwell 
[5] who found that the critical Taylor number, 7’. , was 1707.8 at a. ~ 3.13. 
This boundary value problem is not self-adjoint* and hence can not be expected to 
be represented by a variational principle in the usual manner. However if we let Z = T°”” 
and iaZu = (D*° — a’)v then Eqs. (1) and (2) become 


D? — ayu = iar, (D’ — av = iaZu (3) 


and 
u= Du=v=0 at x= +}. 
Now the boundary value problem defined by Eqs. (3) and (4) is equivalent to the follow- 
ing variational problem: Let 
1/2 


1/2 
ni/ae 


I(u,v) = 3 | {{(D? — aduj’ — (Dv) — a’v’} dx — iaZ uv dx. (5) 
2 J-1/2 


Then among all functions u and v which are continuous and have continuous derivatives 
of the fourth order and second order respectively, and satisfy the boundary conditions 
(4) that pair which makes J stationary necessarily satisfies the differential equations 
(3). To prove this statement we note that the variation in J due to a variation 6 wu in 


u and év in v is 


- 


al/2 


6] = | {(D? — a’)u(D* — a’) bu — iaZv bu} dx 


_ | {Dv D &v + a’v &v + iaZu sv} dx. 
v-—1/2 
If 67 is to vanish for arbitrary variations 6u and év that satisfy the boundary conditions 
(4) it follows immediately after integration by parts that wu and v satisfy the differential 
equations (3). 
To solve this variational problem we first note that from the form of the functional 
I and the boundary conditions (4) that the solution can be split into even and odd func- 
*Briefly, the boundary value problem L(u) = 0 on x; < x < 22, where L is a differential operator, 
is said to be self-adjoint if for any two functions u and v satisfying the specified boundary conditions 
at z; and z, the integral from z; to x2 of vL(u) — uL(v) vanishes. See [8, p. 59] or [9, pp. 42-44] for a 
more detailed discussion of self-adjoint boundary value problems. 
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tions about « = 0. For the even solution we expand wu and v in a complete set of even 
orthogonal functions on — 1/2 < x < 1/2. Appropriate series are 
u(x) = ) a,E,(x), v(x) = = b,E,(2x), (6) 
n=1 n=1 
where E,(a) = 2'* cos (2n — 1)xrx. Notice that the boundary conditions u = v = 0 at 
x + 1/2 are automatically satisfied by these series. The boundary conditions Du = 0 
at « = + 1/2 introduce the constraining condition 
r= >> (—1)"*'(Qn — Da, = 0 (7) 


n=1 


on the a,. Substituting* the series for uw and v in the expression for J gives 


[ = > 14(A2a> — A,b?) — iaZa,b,}, (8) 
n=1 
where A, = (2n — 1)*x” + a’. The a, and b, are determined by requiring that the 


variation of J with respect to the a, and 6, vanish subject to the constraint T = 0. The 
constraining condition is introduced by the use of a Lagrange multiplier. Thus we set 
the partial derivatives of J — wT with respect to a, and b, = 0. Here u is a Lagrange 
multiplier. This gives two simultaneous linear non-homogeneous equations for a, and 6,. 
Solving for a, and substituting in the condition ! = 0 gives the following transcendental 


equation for 7 as a function of a 


(2n — 1)’ yo , 
—{~—* = 0. 9 
> AS —eTr (9) 
Although the series in Eq. (9) only converges like (2n — 1)~*, the convergence can be 


improved quite easily by making use of the fact that the sum of the squares of the 
reciprocals of the odd integers is r’/8. This gives 


> f(2n — 1)°A, ‘oa 
st x {os aT - aia} <0 (10) 


A. -— eT 





and the convergence of this series is like (2n — 1)~*. For a given value of a the roots of 
Eq. (10) can be conveniently found by trial and error. The smallest positive value of T 
is 1708.1 and occurs for a ~ 3.12. The solution for u and v odd can be carried through 
in exactly the same manner. For this case an appropriate set of expansion functions 
would be F(x) = 2'” sin 2nrz. 

It is perhaps worth mentioning that if the Taylor problem is written in the form of 
Eqs. (3) and (4) rather than the original form of Eqs. (1) and (2), approximate solutions 
can be obtained conveniently by the Galerkin method. This method consists of ex- 
panding u and v in sets of complete functions (preferably orthogonal) that satisfy the 
boundary conditions and then requiring the error in the equations for u and v to be 
orthogonal to the expansion functions for u and v. (For a more detailed discussion of the 
Galerkin method, see Chap. 4 of [10]). For the Taylor problem appropriate functions 
are the E(x), mentioned earlier, for v(x) and the set of orthonormal even functions, 


*It is permissible to differentiate these series termwise. This will be discussed in Sec. 4. 
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C,(x)*, satisfying C, and DC, = 0 at x = + } recently tabulated by Harris and Reid 
[11] for u(x). In particular if one term in each series is used we obtain the approximate 
result 7’, = 1728 at a = 3.12. This result can be improved by taking more terms in the 
series; however J/ terms in each of the series for u and v will lead to a determinantal 
equation of order 2M for T as a function of a, and the exact solution will require the 
evaluation of determinant of infinite order. It might also be mentioned that the higher 
characteristic values can be computed easily from Eq. (10). In contrast in the Galerkin 
method, the number of terms required in the approximating series increases as the 
successive characteristic values are computed. 

Thus for the Taylor problem it would appear to be preferable to use the methods 
suggested earlier rather than to use series that satisfy the boundary conditions term- 
wise. In Sec. 5 the use of essentially the Galerkin method to solve exactly the Taylor 
problem as given by Eqs. (1) and (2) by using series which do not satisfy the boundary 
conditions termwise will be discussed. In many problems, however, the methods of 
this paper are not applicable. In particular, this is true for a differential equation with 
variable coefficients and in these cases the Galerkin method may be very appropriate. 
This depends primarily upon the availability of functions satisfying the necessary 
boundary conditions and how closely such functions can be expected to approximate 
the characteristic functions of the problem. (See [12] for a discussion of a particular 
problem). 

3. The inhibition of convection by a magnetic field. The differential equation 
governing the instability of a layer of an electrically conducting fluid heated from below 
in the presence of an external magnetic field parallel to the gravitational field is 

(D’ — a’){(D® — a’)? —Q D*}W = —-akRwW, —4t<2z< i, (11) 
where D = d/dx, a is a dimensionless wave number, Q is a measure of the magnetic 
field strength, and R, called the Rayleigh number, depends upon the temperature 
gradient. The boundary conditions are 


W =0, ((D? —a) —QD’}W =0 (12) 


at x = + 3 and 
DW=0 or DW=0 (13) 


at a rigid or a free surface respectively. (See [3] for a derivation of this boundary value 
problem). The physical problem requires for a given value of Q (real and positive) the 
determination of the minimum real positive value of F for all real positive a for which 
this boundary value problem has a solution. 

The methods of Sec. 2 are applicable to this problem. First let Z’ = R, and define 
U(x) by iaZU(x) = {(D? — a’)? — QD’}W, then Eq. (11) becomes 


{((D’ — a’)’ — Q D*}W = iaZU, (14) 
(D? — a’®)U = iaZW. | 


*The functions C,(xz) are of the form (cosh \mx)/(cosh 44m) — (COS Amz)/(cos 4A»), where the A, 
are the positive roots of tanh 4 + tan 4A = 0. 
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The boundary conditions are 


W=U=0 (15) 
at x = + } and 
1. Two rigid surfaces: 
DW=0 at x= +} (16a) 
2. Two free surfaces: 
D’W =0 at x= +} (16b) 


3. One rigid surface and one free surface: 
DW =0 at x=}, DW=0 at x= —}. (16¢) 


lor the variational problem the appropriate functional is 


KU,W)=43[  {{(D? -— a) Wy + Q(DW) — (DU) = aU} dex 
- (17) 


al/2 
—iaZ | WU adr. 

v-1/2 
For Case 1 it can be shown, as was done in Sec. 2, that a necessary condition for the 
vanishing of the first variation in /(U, W) subject to arbitrary variations in U and W 
that satisfy the boundary conditions (15) and (16a) is that U and W satisfy the differ- 
ential equations (14). This problem is precisely the same as the Taylor problem except 
for the term Q(DW)° in the definition of 7. The appearance of this term does not cause 

any complications and hence this case will not be discussed further. 

For Case 2 the statement of the variational problem is slightly different. The boundary 
conditions D°W 0 at « = + 3} are natural boundary conditions. That is, among all 
functions U and W satisfying the boundary conditions (15) that set which makes J 
stationary will necessarily satisfy the differential equations (14) and the boundary 
conditions D>W = 0 at « = + }. This follows immediately upon noting that 


21/2 


ve ‘(D? — a?)W(D? — a) bW + Q DW D sW — DU DU — aU 8U 
— iaZ(W 6U + U 8W)} dx 


nl /2 


= [D’W D 6W]'2, + |  {(D? — a)"W — Q DW — iaZU} BW ax 


/—1/2 
+ [  ((D* = aU — iaZW) BU az 


after integration by parts and the use of the boundary conditions. The vanishing of 67 
for 6U and 6W arbitrary in — } < x < }. and DéW arbitrary at x = + 3 gives the desired 
result. The computations are particularly simple for this case and will not be discussed. 
Indeed, for a solution that is even about x = 0, one can choose W = a, cos (2n — 1)xz, 
LU) = b, cos (2n — 1)xx. The boundary conditions are then satisfied and the relationship 
between FR, a, and Q is found by substituting in the differential equations and setting 
the determinant of the coefficients of a, and b, equal to zero. 
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For Case 3 the variational principle is the same as in Case 2 except that the con- 
dition DW = 0 at x = 3 must be satisfied by the variation in W; the boundary con- 
dition D?W = Oat x = — #isstill a natural boundary condition. Because of the boundary 


conditions the solution for this case can not be split into even and odd solutions. Conse- 


quently we write 


co 


Wa) = DO {a,£.(a) + eF.(a)}, Ue) = DL {b.EF.(z) +d,F,(a)}, (18) 
= n 1 
where E,(x) = 2” cos (2n — 1)xx and F,,(x) = 2’* sin 2nwx. The boundary conditions 


W = U = Oatx = + } are automatically satisfied by this choice of expansion functions. 


The boundary condition DW = 0 at x = 3 introduces the constraint 
Tr = > (—1)"*"{(2n — 1a, — 2ne,} = 0. (19) 
n=1 


Substituting the series for W and U in the expression for J and carrying out the necessary 


integrations gives J as a function of the a, , b, , c, and d, . The vanishing of the first 


variation in J subject to the constraint (19) requires the vanishing of the partial deriv- 
c, and d, . Here uw. is a Lagrange multiplier. 
c. @. 


n 


atives of J — u,T, with respect to the a, , b, , 
This leads to four simultaneous linear non-homogeneous equations for a, , 6, , 


which can be solved quite easily. Substituting for the a, and c, in the condition I, 0 
gives the following transcendental equation for FR as a function of a and Q 
n?(n*x? + a’) 
Bo — 0. (20) 





> ‘n'e +a’) + nr Q\(n'o +a) —aR 


Although this series only converges like n~ its convergence can be improved by adding 
and subtracting 7°/6 = > 03 n~ 
Numerical computations were carried out for two cases. In the case Q = 0 ten terms 


of the series were used to compute R, = 1100.6 at a ~ 2.86. This case was also solved 
1100.65. For 


by standard techniques by Pellew and Southwell [5] who found R, = 
Q = 50, a = 3.45 fourteen terms in the series were used to determine R = 2217.4. It is 


estimated the error is less than .5%. For this case Chandrasekhar [3] gives an upper- 
bound of 2217.6 for R*. 

4. The Taylor problem with an axial magnetic field. 
discussed in Sec. 2 with the additional complicating feature that the fluid is an electrical 
conductor and there is an axial magnetic field. In this case, the stability problem requires 


Consider the Taylor problem 


the solution of 
(21) 


bol 
tol 


{(D® — a’)? + Qa’}*v = —Ta*(D’ — a’), —-i<z< 


with the boundary conditions 
2 a ) 
Dv = (D° — a )v = 0,7 | (22) 
{(D? — a’)? + Qa’}\v = 0, D{(D® — a’)? + Qa’}v = OJ 


at z = + 3. Here Q is a measure of the magnetic field strength, D = d/dz and a and T 


*In Table 3 of [3] the values of Q should be divided by four and the values of a should be divided 


by two as pointed out by Chandrasekhar on p. 1187 of [4]. 
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have the same meaning as in Sec. 2. (See [2] for a derivation of Eqs. (21) and (22).) 
Again, for an assigned value of Q, we are interested in determining the minimum positive 
real value of 7’ as a function of a (positive) for which this boundary value problem has a 
solution. 

To formulate this problem as a variational problem, we again split the differential 
equation into two lower order equations. Also in this case it is convenient, though 
not necessary, to change the interval of integration. If we let x = x(z + 3), b = a/z, 
P = Q/r, Z = T/r* and ibZu = {(D’ — b’)? + Pb’}v, where D now denotes d/dz, 


Eqs. (21) and (22) may be written as 
{(D? — b?)? + Pb?}u = ibZ(D? — bv (23) 
{(D? — b’)? + Pb?}v = ibZu, O<a<fr 


and 
u= Du = Dv = (D’ — wv = 0 (24) 


at x = O and z. In contrast to the problems treated in Secs. 2 and 3 the above boundary 
value problem, even after this splitting, is not self-adjoint and hence does not have a 
variational principle of the type derived previously. 

In order to use the methods of the calculus of variations, it is necessary to introduce 
the adjoint boundary value problem. To do this, we write Eqs. (23) as L(u) where uw is 
the vector (u, v) and L is the matrix differential operator defined by Eqs. (23). We then 
form the integral of UL(u) from 0 to x where U = (U, V). With the use of integrations 
by parts this integral is then written as an integral from 0 to z of uL,(U). The adjoint 
differential equations are L,(U) = 0. The boundary conditions are determined by requir- 
ing that the terms resulting from the integrations by parts vanish. (See p. 149 of [9]). 
For this particular problem the adjoint problem is 


{((D? — b)° + PUIU = ibZV (25) 
{(D® — b’)? + Pb’}V = ibZ(D? — b*)U, O0O<a<fr 
and 
U = DU = V = D(D’ — W)V =0 (26) 
at x = 0 and z. Notice that if U and V are associated with u and v the boundary con- 


ditions on V are different and also the role of the operator (D’ — b’) is reversed. 
The original boundary value problem and the adjoint boundary value problem 


have the following variational principle. Let 


I(u,v, U, V) = [ {M(u)M(U) + M@)M(V) + Pb’(uU + vV) (27) 
+ ibZ(Dv DU + b*v0U — uV)} dz, 
where M(u) is (D’? — b’)u. Among all functions u, v, U, and V that are continuous and 


have continuous derivatives through the fourth order and satisfy the boundary con- 
ditions. 
u= Du=Dv=0, U=DU=V=0 (28) 


at x = 0 and z that set which make the first variation in J vanish necessarily satisfy 
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the differential equations (23) and (25) and the natural boundary conditions 
(D? — b’)v = 0 and D(D? — b’)V = O at the end points. The proof of this variational 
principle is straightforward and will not be given. 

The solution of this variational problem can be split into non-combining even and 
odd (about +/2) solutions. For the even solution appropriate orthogonal series are 


U= A+ bs a, COS NZ, U=b+ b ¥ b,, cos ~“ (29) 
v= + den COs NZ, V=da+ > d, cos nz 


where >, stands for summation over the even integers. Although the series for u, v, and 
U can be differentiated termwise twice, this is not true for the series for V. To see this 
we note that since V is even, DV is odd and D’V is even; hence appropriate series for 
these functions are DV = > dé sin nx, and D’V = di’ + >> d?’ cos nx. To determine 
the relationship between the d, , d’ and d/’ we integrate by parts the following expressions 
for di and d’ 


rT 


r 9 
ad, = 2/ DV sin nz dz = = Ven nzjo ——n / V cos nx dx (30a) 
TFJo Tv Tv ry) 
= —ndn, n= 2,4,6- 
and 
rt ] , 2 l ri" 1 
Pa - | D*V dz = =[DV]< = 36 (30b) 
T Jo T 
say, and 
2 ; 217 o : _2 7 1 * 2 . 7 2} 
d,’ = - [ D°V cos nz dz = 2 [DV cos nz]o + =n | DV sin nx dx (30c) 


=B+nd,=8-—n'd,, n=2,4,6--:. 


Similar computations show that the series for u, v, and U may be differentiated term- 
wise twice. The idea of expanding a function and its derivatives in series of complete 
functions and then integrating backwards to obtain the relationship between coefficients 
is not new. For instance it was used by 8. Goldstein in 1936 [13]. 

Once appropriate orthogonal series for u, v, U and V and their derivatives have 
been obtained the solution of the variational problem, while slightly more complicated 
than in the previous cases, is straight forward and hence will only be briefly described. 
The boundary conditions Du = Dv = DU = 0 at x = 0 and « are automatically satis- 
fied by the choice of the expansion functions. The boundary conditions u = U = V = 0 
at x = 0 and = introduce the constraining conditions. 


Ir=a+)>4,=0, Kw=b+>b=0, T;=a+ Dd =0. (31) 


Substituting the series for u,v, U and V in the expression for J and making use of relation- 
ships between d/ , d’’ and d, gives J as a function of the a, , b, , c, , d, and 8. The vanish- 
ing of the first variation in J subject to the constraining conditions (31) requires that the 
partial derivatives of J — u3;0; — ul, — usI's with respect to the a, , b, , c, , d, be zero. 
Here py; , u, , and ws are Lagrange multipliers. This gives four simultaneous linear non- 
homogeneous equations for the a, , b, , ¢, , d, which can be solved to give these co- 
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efficients as linear functions of ws , 4, , us and 6. The fact that we do not obtain infinitely 
many equations is a result of using series with the appropriate orthogonality properties. 
Substituting for the a, , b, , and d, in the constraining conditions (31) and the additional 
necessary condition for the vanishing of 67, 0(I — w3T3 — usT', — usl's)/08 = 0, gives 
four simultaneous linear homogeneous equations for us , u, , Hs and 8. The vanishing of 
the determinant of the coefficients gives a four by four determinantal equation for Z 
as a function of b and P. This determinant reduces immediately to a two by two determi- 
nant which can be written as 


ims =. ey O = 0, (32) 
where 
- _ (bd? + P) Br + Pb° ss B, 
es +2> os » BeztrIze, 
4 2 2 2 
Sa b°(b> + P) 42 y B,(B, + Pb) 
Ao A, 


and B, = n? + b?, A, = b* {(b? + P)’? — 2°}, 4, = (B2 + Pb’) — b°Z’B, . Although 
the series in X, only converges like n~* its convergence can be improved by adding and 
subtracting 7’/24 = >> n~*. Now all the series converge like n~* where summation is 
over the even integers, hence only a few terms are needed to evaluate these 
series accurately. 

In the case P = 0 it can be shown that this result reduces to the solution of the 
Taylor problem discussed in Sec. 2. Of course in this case these series are not as convenient 
as those used in See. 2. For the case Q = 100, a = 3.35 a sample computation was carried 
out using six terms in the series. The smallest positive root was found to be 
T = 1.759 X 10*. It is estimated this result is correct to less than 1 %. Chandrasekhar 
[3] gives JT = 1.757 X 10* for this case. 

5. Direct series solution. Although we have used variational techniques to solve 
the problems discussed in the previous sections, this is not necessary. Indeed in some 
cases it may be more convenient to use the direct method illustrated below for the 
Taylor problem that was discussed in Sec. 2. 

First introduce the new variables z = x (x + 4), b = a/x, and Z’ = T/z* so that 
the boundary value problem defined by Eqs. (1) and (2) becomes 


(-iva 0H, O<s<e (33) 
v=(D? —bwv= D('D? —VUwv=0 at z=O0 and fz, (34) 


where D now stands for d/dz. For the even solution (about z = 2/2) of Eqs. (33) and (34), 
we expand v and its derivatives in the following series. 


v= > a,sin nz Dv = > e, sin nz 
Dv = >> b, cosne Dv = > f, cosnz 
Dv = Doe, sin nz Dv = > g, sin nz 
Dv = > d, cos nz 


where >. now stands for a summation over the odd integers. This notation will be used 
throughout this section. Substituting these series in the differential equation (33), 
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multiplying by sin mz, m = 1, 3, 5, --+ and integrating from 0 to x gives 
In — 3b’e, + 3b‘c, —_ b°a,, = — b*Z'a, : n= 1.3,5->-. (35) 


This is similar to the Galerkin method in that the error in the differential equation is 
required to be orthogonal to the expansion functions; however it should be noted that 
the series used here do not satisfy the boundary conditions termwise. Indeed although 
the first two boundary conditions are automatically satisfied the boundary condition 
D(D* — b’)v = 0 at z = O and = requires 


> (d, — n’b,) = 0. (36) 
The relationships between the coefficients in the series for v and its derivatives are 
found by integrating backwards as was done in the previous section. Thus, }, = na, , 
¢, = —na,,d, = — n'a,,e, = + n'a, , and 
of Ss 2 , Aeron ee a ; 
i=" Dv cos nz dz = —4§[D'v cos nz]p +7 Dv sin nz dz? = —— D'v(0) + ne, 
T Jo T 0 vs 
orf, =y + ne, =y + na, say, and g, = — ny — na, . In deriving these relationships, 
use has been made of the boundary conditions that v satisfies. Substituting for g, , ¢, , 
and c, in Eq. (35) and solving for a, gives a, = — ny/{(n*> + b°)* — b°Z*}. Also Eq. 


(36), expressed in terms of the a, is, >> n(n? + b*)a, = 0. Substituting for the a, gives 
the following transcendental equation for Z as a function of b, 
n*(n? + b?) _ — 
2a + oy — 07 =O 30 
This equation is identical with Eq. (9) derived in Sec. 1 using a variational procedure. 

It should be noted that a, ~ n~° and hence the solution of the characteristic value 
problem can be differentiated termwise four times. This result is in contrast to the 
remark made by Jeffreys that the trigonometric series obtained by differentiating more 
than twice are divergent [18]. Further a simple computation shows that f, ~ n ~~ and 
hence the series for D*v and D°v are convergent. 

6. Remarks. In the previous sections, general methods for solving exactly constant 
coefficient characteristic values problems have been illustrated. The methods depend 
primarily on the expansion of the unknown function and its derivatives in sets of appropri- 
ate orthogonal functions. It is not clear that a general conclusion can be drawn about 
the advantages or disadvantages of the variational method as used in Secs. 2, 3, and 4 
compared to the direct method illustrated in Sec. 5. However it should be mentioned 
that if the direct method is used it is necessary to enforce all the boundary conditions 
even if some of these boundary conditions are natural boundary conditions for the 
rariational problem. 

A partial list of other hydrodynamic and hydromagnetic stability characteristic 
value problems that can be treated easily and exactly by the methods outlined above are 
the inhibition of convection by a magnetic field for a layer of fluid heated from below 
when the magnetic field is not parallel to the gravitational field [4], the stability of a 
viscous flow between two concentric cylinders rotating in the same direction in the 
presence of a circular magnetic field when the gap between the cylinders is small [14], 
and the stability of liquid helium II between concentric rotating cylinders in the case 
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that the cylinders are rotating in the same direction and the gap between the cylinders 


is small [15, 16]. 


These methods have also been used to treat problems in the vibrations of beams and 


several two dimensional problems in elasticity. In particular, the reader is referred to 


references [17] and [19]. 
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BOOK REVIEWS 


(Continued from p. 874) 


Professor Stiefel’s introduction illustrates that the self-adjointness of a continuous problem is most 
naturally preserved in its discrete analogue, not by approximating a partial differential equation and 
its boundary conditions, but by approximating and minimizing the underlying energy integral. Rutis- 
hauser’s essay summarizes the theory of the gradient methods to be tested, and Ginsburg reports on 
the results of his experiments. Engeli deals compactly with both the theory and the numerical conse- 
quences of overrelaxation. 

The verdict, set down in the last chapter, is generally speaking for overrelaxation on Dirichlet- 
type problems and elimination otherwise; for eigenvalue problems, a champion among the gradient 
methods is designated, but the authors make it plain that their investigation did not include some 
other competitive methods. 

The clarity is unimpaired by minor errors in spelling and usage, but a misplaced formula on pp. 
52-53 is more serious, and the notations cosh and Ch both appear. 

This example of a valuable service which can be performed ‘‘collectively’’ by a computation labo- 


ratory unquestionably deserves to be followed. 
W. GILBERT STRANG 


Probability and statistics. Edited by Ulf Grenander. The Harald Cramér Volume. Almqvist 
& Wiksell, Stockholm, and John Wiley & Sons, New York, 1959. 434 pp. $12.50. 


This book contains nineteen papers dedicated to Harald Cramér in honor of his 65th birthday. 
Most of the contributions are survey papers of recent developments and collectively they cover a wide 
range of selected topics in mathematical statistics and probability. The contributors are T. W. Ander- 
son, M.S. Bartlett, J. L. Doob, G. Elfving, W. Feller, E. Fix, J. L. Hodges, E. L. Lehmann, U. Gren- 
ander, M. Kac, D. G. Kendall, P. Lévy, P. Masani, N. Wiener, J. Neyman, H. Robbins, M. Rosenblatt, 
C. O. Segerdahl, T. W. Tukey, S. S. Wilks, and H. Wold. The quality of the papers is as one would 


expect from such a distinguished group of authors and is a fitting tribute to Cramér. 
G. F. NEWELL 


Matrices. By William Vann Parker and James Clifton Eaves. The Ronald Press Co., 
New York, 1960. viii + 195 pp. $7.50. 


This is a textbook on elementary matrix algebra for introductory college courses, probably most 
suitable for students majoring in applications of mathematics. It is clearly written and produced, and 


includes many worked and unworked examples. 
W. FREIBERGER 


Introduction to mathematical statistics. By Robert V. Hogg and Allen T. Craig. The 
Macmillan Co., New York, 1959. ix + 245 pp. $6.75. 


This is one of several new text books on mathematical statistics to appear this year apparently as 
a consequence of the trend toward inclusion of this subject in the mathematics curriculum of most 
American universities at the advanced undergraduate level. The book covers an exceptionally wide 
range of topics including such things as sufficient statistics, most powerful tests and other topics not 
always included in an elementary book. To accomplish this, however, the author must sacrifice mathe- 
matical elegance. Some theorems from other branches of mathematics are quoted without proof, and 
the proofs of some key theorems in probability theory are incomplete. Some people will also object to 
the use of a probability density function for discrete and continuous random variables instead of the 


(Continued on p. 394) 
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SOME COMMENTS ON NARROW BAND-PASS FILTERS* 


BY 
,M. ROSENBLATT 


Brown University 


1. Introduction. It appears to be part of the engineering folklore that a narrow 
band-pass filter applied to a stationary random input yields an output that is approxi- 
mately normally distributed. Of course, such a conjectured result could not be true in 
absolute generality. At the very least, one ought to require ergodicity of the random 
process being filtered. However, we should like to sketch out a domain within which 
such a result would in fact hold and indicate roughly boundaries outside of which such 
normality would not be expected. 

2. Preliminary discussion. Consider x(t), —7 < t < ©, a strictly stationary 
separable (see Doob [1] p. 52) random process. The separability is assumed so that 
meaningful statements can be made about sample functions (or realizations) of the 
process. Later on, a mixing condition and various moment conditions will be introduced 
and assumed to hold for x(t). Our interest will be centered about weighted averages 


aT 


J 


of a certain form. Though we shall actually deal with a more general context, the weight 
function w(t) can be thought of as the weighting introduced by the filter. The assump- 
tions on w,(t) will be such that it is amenable to a generalized Fourier analysis (see 


[2] p. 233). Set 


w7(t)a(t) dt (2.1) 


i 


W=Wwer)=[ wi(o de. (2.2) 


An explicit indication of the dependence of W on 7’ will often be deleted in the rest of 
the paper 
The conditions on w(t) will be called assumptions (A). 

Assumptions (A 

1 W(T) > e©asT> & (2.3) 
2. The functions wr(t) are slowly increasing in that 


a | w(t) >dt = o(W) as TT @ (2.4) 
“A(T 
for any sequence of subsets A(T) of [0, T| with the Lebesque measure of A(t) m(A(T)) = 
o(T), uniformly in m(A(T))/T, as T — © and 


b) w(t) = O(W'”) uniformly in t (2.5) 


*Received March 14, 1960. The work presented in this paper was carried out on behalf of the Bell 
Telephone Laboratories. The author gratefully acknowledges their support and encouragement. 
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as T — 
a T-—\h\ 


3. lim W | wr(t + | h |)wr(t) dt = p(h) (2.6) 
T-2 70 

exists for every h and is continuous tn h. 

The limit function p(A) can be seen to be a non-negative definite function and therefore 

by Bochner’s theorem ({4] p. 207) it has a representation 


p(h) = | e'™ dM(a) (2.7) 
with M(A) a non-decreasing bounded function. The spectral mass of M(A) is sym- 
metrically located about zero since p(h) = p(— h). When w(t) corresponds to a band- 
pass filter, conditions 1 and 2 of (A) will usually be obviously satisfied since the functions 
w(t) are typically uniformly bounded in 7 with W(7') of the same order of magnitude 
as 7’. If one is narrow band-pass filtering about yu, 1/7(A) will increase only at u and — ux. 
The conditions on w7(t) have been taken somewhat more general than required in band- 
pass filtering so as to cover situations that commonly occur in discussions of regression 
problems as they arise in time series analysis ((2], p. 233). 

The process x(t) will be assumed to satisfy a strong mixing condition. Let @, be the 
Borel field of events generated by the random variables x(u), u < t, and $, the Borel 
field of events generated by the random variables x(u), u > 7. Thus ®, and §, represent 
the information given by knowledge of the process before t and after 7 respectively. 

Assumption (B). The process x(t) satisfies the strong mixing condition (B) if there is 
some positive function g(u) defined forO0 < u < @ withg(u) | Oasu— © such that for 
any pair of events Be2G, ,? <3, .t < =. 


| P(B CO F) — P(B)P(F) | < g(r — 2). (2.8) 


The mixing condition (B) might be called a uniform mixing condition. It is appreciably 
stronger than ergodicity and somewhat stronger than what one ordinarily calls mixing 
((3] p. 36). We shall later see by virtue of an example that ergodicity alone is not enough 
to guarantee that the output of a narrow band-pass filter be normal. 

As in the case of any central limit theorem, certain moment conditions are necessary. 
The conditions to be specified are certainly not the best that could be obtained but they 
are quite natural in the context we are talking about. Assume that fourth moments 
Ex(t)* exist and that the process x(t) is continuous in the mean of fourth order, that is, 
— 0 (2.9) 


4 


E | x(t) — 2x(z) 


Because of the stationary of x(t), the mixed fourth order moment 


E[x(t,)x(t.)x(ts)a(t,)] = P(t. — th , ts; — th ty — th) (2.10) 


as t — r. It is convenient to take the constant first moment F x(t) as identically zero. 


depends only on the differences t, — t, , ts — t, , ts — t, . Let R(t) denote the covariance 
R(t, — t,) = Elx(t,)x(t.)}. (2.11) 
We introduce the fourth order cumulant function 
Ot ang t; » ts a t; ’ ts — t,) = P(t, — ty » ls — t »t aa t,) (2.12) 
—_ Po(t — ty ; ts aed ty ; te — t,), 
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where P, is what P would be in the case of a normal process, namely 
Po(te — ty ; t, ate t, ; ts — t;) = R(t. — i) R(t, — t;) a R(t; — t, R(t, = t) 
+ R(t, — t)R(t; — te) 


(see Magness’ paper [5]). Notice that if R(¢) is absolutely integrable, the spectral density 
f(A) of the process exists and is continuous. 

Assumption (C). R(t) and Q(t, , te , ts) are absolutely integrable over one and three 
dimensional space respectively. The spectral density f(A) ts assumed to be positive every- 
where. If one were only considering narrow band-pass filtering about u, the positivity of 
{(\) would only be required at yu. It should also be noted that the proof of asymptotic 
normality to be given does not use the stationarity of x(¢) fully, only stationarity up to 
moments of fourth order. 

3. The basic theorem. The basic scheme of the proof of the theorem follows that 
given in [6] and is a formalization of a heuristic notion of A. Markoff. The derivation 
of the result will be given in detail and in several stages, even though this will duplicate 
much of the discussion of [6]. 

Theorem. Under assumptions (A), (B), and (C) 


- 
ie [ w7(t)a(t) dt (3.1) 
is asymptotically normally distributed with mean zero and variance 
Qn [ f(s) dM(). (3.2) 
Let 
7 
S(T) = [ wy(t)x(t) dt (3.3) 
and set 
a (itl) p(T) +ia(T) 
U(T) = | w r(t)a(t) dt, (3.4) 
Ji(p(T)+e(T)] 
o(i+1) (p(T) +e(T)) 
VAT) = | w (t)2x() dl, (3.5) 
/ (+1) p(T) +ie(T) 


j = 0,1, ---,k — 1. Here k[p(T) + q(T)] = T and the numbers k = k(T), p(T), q(T) 
will be chosen so that k, p,q — © and q(T)/p(T) — O0as T — ~. In effect, the interval 
[0, 7] is being divided into an alternating succession of big blocks (length p(7)) and 
small blocks (length q(7’)). The U;,’s are the large block integrals and the V;’s the small 
block integrals. We first show that the contribution of the small block integrals is negli- 
gible. Now 
k 2 | p 

E|> VW”?! =E| | w (a(t) dtw'? 
9=1 J A(T) 


|2 
| 
| 

| | 


< [ I r(u) | [ | w(t + u)we(t) | dt du W(T)" (3.6) 
ao JA(T) 


lA 


/ | r(u) | if. | w(t + u) |? dt [.. | w(t) |’ a) “war, 
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where 
A(T) = U {t| wi) + G - Nal) < t < iT) + A(T). 
j=1 
But (3.6) must approach zero because of the absolute integrability of r(u) and part 2a 
of assumption (A). But this means that 


k 
> Vw"? +0 (: 
i 


~I 


ve 


in probability as T ; 
Our object now is to show that 


> UW”? (3.8) 


. > CO 


is asymptotically normally distributed with mean zero and variance (3.2) as 7’ 
The theorem on asymptotic normality of S(7)W~'” would then follow immediately 
for (3.7) is asymptotically negligible. Let us introduce for this purpose the distribution 


functions 


G;.r(z) = P{U;W” < z} (3.9) 
and the events 
A(j, 77, m; , 6) = {m; 6 < U;,W'” < (m; + 1) 3}. (3.10) 
Now 
k k 
(a A(j, T, m; , s)) < o> ta" =< x) 
(my teetemetkhyicz \ je] j=] 
(3.11) 
< be (A Aj, T, m; ; i)). 
Notice that 
P( max |U;W'? | > \n) <e, (3.12) 
j=1,°**,k 
where 7, = (kC/e)'” with C a constant. This follows from the fact that 
| oe 2 | k j2 
E| max |UW'’|| <E) D | Uw’ || 
| fmt,eoe,k ie (3.13) 
k 2 
< weary k'” U; P) < LC 
and an application of the Tchebycheff inequality. 
Further 
| . , 
| be (1) A(j, T, m; , i)) _ a [] P(AG, 7, m; , 4)) 
5< 1 (miters+me)dsz 1 | 2 
3.14) 


< u( 724) g(a) + 2c. 
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For the probability contributed by all the sets  * A(j, T, m; , 6) for which max 
U;,w'”"| > 7 is by (3.12) at most e. Consider now all the sets (\;-. AG, T, m; , 8) 
for which max |U;W~'”’| < 7, . By repeated application of the mixing condition (B) it 
is clear that 

| k k 

P(A AG, 7, m8) — TI PAG, 7, m8] <koQM). 8.18) 


| i=1 


Since there are (27,/6)* sets of this form, the desired inequality (3.14) is obtained. 

Inequality (3.14) will later be applied to show that the sum of the U;W~’” has the 
same asymptotic distribution as the sum of independent random variables with the 
same marginal distributions, as long as k(7'), q(T), p(T) are appropriately chosen. 
However, let us see what the asymptotic distribution of the sum of such independent 
random variables would be. The distribution function of these k independent random 
variables is 


Gi,7r* +++ * Gy, r(2), (3.16) 


the convolution of G,,r(x), «++ , Ge, r(%). 
Lemma. The distribution function (3.16) is asymptotically normally distributed with 


mean zero and variance (3.2). 
Now 


k o 
DS E| Uw |? — | r(u) p(u) du 


(3.17) 


a 


= Se [ f(s) dMQ) 


as T — @ by parts 2a and 3 of assumption (A). Further (3.17) is positive since f(A) 
is positive everywhere. By Liapounov’s form of the central limit theorem [4], if 


2 


k k 
DS E\|UW’ | = o( |U,w’” ’) (3.18) 
j=l i=l 


the expression (3.16) is asymptotically normal with mean zero and variance (3.2). But 


(i+1) p(T) +ia(T) 


E| Uw” |* = WT)” il W r(t,) Ww r( te) w r( ts) w r( ty) 


i(p(7T)+¢e(T)) 
[R(t — t,)R(ts a t,) + R(t, — ts) R(t nat t,) 
+ Rit, — t)R(tz — tb) +QOb—-th,t—-h,u- t,)| dt, dt, dts dt, . 


The sum over j of the first three terms on the right hand side of equality (3.19) is 


k 
3 DE’ | U,w-” |’. (3.20) 
i=1 


(3.19) 


By part 2a of assumption (A), Z |U;W~'”’/? approaches zero uniformly in j. This coupled 
with (3.17) implies that expression (3.20) approaches zero as JT’ — ~. Consider now the 
last term on the right hand side of (3.19). By part 2 of assumption (A) and the absolute 
integrability of Q, the sum over j of the last term on the right of (3.19) tends to zero as 
7 — o., Liapounov’s condition for the central limit theorem is therefore satisfied and 


the Lemma is established. 
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Notice that 


k 
Girt *G,7e@ —kd < = T] P(AG, 7, m; , 8)) 
i7=1 
(my, +***+mytk) bez 


< : 2 I] P(AG, 7, m; , 8) S Gir * +++ * Gere + k 8). (3.21) 


. +m,)i<z 


The Lemma coupled with (3.7), (3.11), (3.14) and (3.21) implies the desired theorem 
if 6(7), k(T), p(T), q(T) can be chosen so that 


K(T)[p(7') + 9(T)] = T 
kK(T), p(T), oT) > © 

onsite (3.22) 
k(T) 6(T) +0 





(2) g(q(T')) — 6. 


The condition k(7)6(T) — 0 is easily satisfied if we set 6 = k-*. The difficult condition 


to satisfy is the last one. Now 


k(2r,/6)* < k™* D', (3.23) 
where D = 2(C/e)'’~. Given the existence of a function g satisfying assumption (B), 
it can always be taken so that 
gu) > (u+ 1)". (3.24) 
If k is chosen so that 
k < [—log g(q(T))]}'” (3.25) 


the last of the conditions (3.22) is satisfied. Keeping these remarks in mind, it is clear 
that if one takes g(7’) = T'”’ for large 7’, all the conditions (3.22) can be satisfied. The 
proof of asymptotic normality is complete. 

4. Some examples. As noted in the introduction, the general expectation is that 
the output of a narrow band-pass filter ought to be approximately normally distributed. 
The theorem obtained in Sec. 3 indicates a fairly broad domain within which this 
expectation is realized. Our object here is to present a few examples that would roughly 
describe boundaries outside of which approximate normality would not be valid . 

It is clear that the filtered process generally ought to be ergodic at least if the output 
of the filter is to be approximately normal. We present a simple example to show that 
ergodicity alone is not enough. Let f(t) be a bounded continuous function with period 


one and 


[ f(t) cos 2xt dt ¥ 0. (4.1) 


#0 


al 
| f(t) dt = 0, 


Let the process x(t) be given by 
a(t) = f(t+ a), (4.2) 
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where a@ is a random variable uniformly distributed on [0, 1]. The process x(t) is ergodic 


with mean value 


1 
Ex(t) = [ f(t +a) da = 0. (4.3) 
Consider the narrow band-pass filter about 2x 
o7 al T 
S(T) = | a(t) cos 2xi dt = [T] | f(t + a) cos 2rt dt + x(t) cos 2rt dt. (4.4) 
J0 Jo J{(T) 


Here [7] denotes the largest integer less than or equal to 7’. Now 


~T 


lim i | 


T-+@ 


1 
x(t) cos 2rt dt = [ f(t + a) cos 2rt dt (4.5) 


which is non-normal. However, it should be noted that this process x(t) is not mixing 
since it has a jump in its spectrum at 27. 

The second example indicates that if all assumptions of the theorem except for the 
positivity of the spectral density at the frequency at which the filter is set up are satisfied, 
the conclusion of the Theorem may not hold. Let &(¢) be a homogeneous differential 


process with non-normal increments. Take 
a(t) = &(t) — 2&(t — 1) + &(t — 2). (4.6) 


Consider the narrow band-pass filter 


aT 


i x(t) dt (4.7) 


about zero. Now 


id 0 


[ x) dt = [ (E(t) — £(t — 1)] at -f (E(t) — &(t — 1)] dt (4.8) 


-1 


so that the output is not normal. Assumptions (A), (B), (C) are satisfied except in that 
f(A) has a high order contact with zero at A = 0. 

5. Acknowledgment. The problem discussed in these comments was suggested 
by S. P. Lloyd. Much of what is relevant in the comments is based on stimulating dis- 
cussions with S. P. Lloyd and D. Slepian. 
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BOOK REVIEWS 
(Continued fro p. 886) 


distribution function. If one takes the attitude that one should try to cover as many topics as possible 
and not avoid certain subjects just because a good mathematical treatment of them is either too difficult 


or too long, then this would certainly serve as a satisfactory introduction. 
G. F. NEwELL 


Modern probability theory and its applications. By Emanuel Parzen. John Wiley «& 
Sons, Inc., New York and London, 1960. xv + 464 pp. $10.75. 


This book is intended as a text for an undergraduate course in probability theory. It is novel and 
may cause controversy but is sufficiently well written that it is not likely to be ignored. It is perhaps 
the most significant introductory book written on this subject since Feller’s famous book in this same 
Wiley series. 

It is claimed that about three quarters of the book can be covered in 39 lectures to undergraduates 
with one year of calculus. Since the author has undoubtedly done so one must believe that it can be 
done but certainly the average student at this level must find this pace very strenuous. The book seems 
more suited for advanced undergraduates or first year graduate students. 

The style of writing through about the first half of the book reminds one of Feller’s book and reflects 
the author’s comment that “probability theory is a subject with great charm and intrinsic interest of 
its own.” The routine material on combinatorial problems, simple types of distributions, etc. is de- 
scribed as if it were alive and the author makes copious use of examples the answers to which are con- 
trary to intuition. About a third of the way through the book, probabilities on real number sets are 
introduced and the second half of the book deals mainly with continuous random variables including a 
discussion of limit theorems for sums of independent random variables. 

Since ‘modern probability theory” is based upon measure theory and involves rather sophisticated 
theorems, it is not possible to include in an undergraduate text a treatment that is both rigorous and 
reasonably complete. Here the author uses modern terminology and at least sketches the line of reason- 
ing behind the key theorems of probability theory with frequent use of such phrases as “in more ad- 
vanced treatments one can prove—.”’ Although the intuitive meaning of such things as Borel sets and 
random variables are explained, precise definitions are frequently lacking. Thus the book is a compro- 
mise between the treatise that gives detailed proofs of key theorems and the more usual elementary 
text that just avoids them. 

The book contains numerous references including the most recent and in most cases the best avail- 
able books or papers (many of which, however, would be too advanced for an undergraduate to under- 
stand), many exercises and problems, and many interesting historical sketches. It establishes a style 
of undergraduate course that many teachers will wish to follow and perhaps this first edition will be 


the prelude to still better ones. 
Gorpon F, NeEwE.u 


Les Transformations de Mellin et de Hankel. By 8S. Colombo. Centre National de la 
Recherche Scientifique, Paris, France, 1959. 99 pp. $2.05. 


This short book is directed to physicists wishing to learn something of the integral transforms as- 
sociated with the names of Mellin and Hankel. The impression it gives is that it is designed for French- 
speaking scientists, since it does not seem to contain anything not already available in books (in 
English) cited by the author in the bibliography. 

The first chapter is an introduction to the study of integral transforms in general, reference being 
made to the Fourier and Laplace transforms whose principal properties are listed and tables of which 
are provided. Next, the Mellin transforms are introduced, their main properties discussed and their 


(Continued on p. 412) 
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—NOTES— 


A SIMPLIFIED GEOMETRIC PROOF OF THE 
RECIPROCAL STRESS THEOREM* 


By B. PAUL (Brown University) 


The well-known reciprocal stress theorem [1, 2]'*? states that the scalar product 
of the stress vector on a given plane with the unit normal to a second plane equals the 
scalar product of the stress vector on the second plane with the unit normal to the 
first plane. Starting from this single theorem, which Biezeno and Grammel [1] call 
“the essential core of all stress analysis” it is a simple matter to deduce the symmetry 
of the stress tensor, to derive the equations of stress transformation, and to find many 
other useful results. 

Biezeno and Grammel prove the theorem geometrically by projecting the faces 
of a stressed tetrahedron and their corresponding stress vectors orthogonally onto 
an auxiliary plane. This is essentially a proof by a technique of descriptive geometry. 
The importance of the theorem may perhaps justify the small amount of space required 
here to present a simpler modified version of this proof which obviates the need for 
projections onto an auxiliary plane. 

For simplicity consider a homogeneous state of stress where S, and S, , respectively, 
represent the stress vectors on planes II, and II, which are respectively normal to unit 
vectors n, and n, . Figure 1 shows two congruent triangles ABC and ABD constructed. 
on the line of intersection of II, and I, . The four points ABCD enclose a tetrahedron. 
Let E, F, G, H represent the centroids of the four faces as shown in Fig. 1. Then it 
follows that the plane of EGH is parallel to II, and the plane of FGH is parallel to II,. 





Fic. 1. 


*Received Oct. 19, 1959 
1C. B. Biezeno and R. Grammel, ‘Engineering Dynamics, Vol. I,” p. 3, Blackie, 1955. 
21. S. Sokolnikoff, “Mathematical Theory of Elasticity,” 2nd Ed., p. 43, McGraw-Hill, 1956. 
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The resultant force vectors on faces BCD and ACD act at the centroids H and G of 
these faces and therefore contribute nothing to the sum of moments about axis GH. 
Denoting the mid-point of GH by Q, decompose the stress vector S, acting at E into 
three components that are acting along QZ, parallel to QH, and normal to the plane 
EGH (i.e., parallel to n,), respectively. Of these only the last component, which has 
the value S,-n, , produces a moment about axis GH. Similarly, only the component 
of S, , acting at F, in the direction of n, , produces a moment about axis GH. The re- 
sultant moment is given by 


(S, nm.) A(EZQ) — (S.-n,) A(FQ) = 0. 


In the above equation A is the common area of triangles ABC and ABD. Since the 
moment arms EQ and FQ are equal by construction, it follows that 


S,-n. = S.-n, . 


UNSTEADY MOTION OF AN INFINITE LIQUID DUE TO 
THE UNIFORM ROTATION OF A SPHERE, r = a* 


By C. D. GHILDYAL (Lucknow University) 


Introduction. The problem of steady motion of a viscous liquid due to the slow 
rotation of a sphere has already been discussed by various workers in the field of hydro- 
dynamics [1]. Here we propose to discuss the unsteady flow of a liquid initially at rest 
due to the uniform rotation, 2, of a sphere about the axis of z, under an external force 
zw acting per unit mass of the liquid along the axis of rotation and the pressure at any 


point of the liquid is given by 
p= il prw dr + aconstant, (A) 


where w = w(r, t), 7° = x? + y’ + 2’, and velocity components at any point (z, y, 2) 
of the liquid are assumed to be u = — wy,v = wrandw = 0. 
The result obtained is valid also for slow rotation of the sphere. In that case the 
external force does not exist and the pressure remains constant throughout the liquid. 
The equations of motion of a viscous homogeneous incompressible liquid in this 


dw (Zs 4 2) 
en ee 4 SZ) 


case reduce to 


ot or r Or (B) 
Op _ 2 
ar prw . 


A general solution of (B) is obtained by the method of Laplace transform under the 


boundary conditions 
(i) w(a,t) = Q, 


(ii) w(r, 0) = 0, (C) 
(iii) limo(r, t) = 0. 


r+ © 


*Received August 26, 1959; revised manuscript received January 20, 1960. 
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T and «(r, T) = Wr, T)/r°”, in (B) and (C), we get 
ow 9W (1.1) 


1. Solution. Putting vt = 
aw _ aw 
oT =—s sar” 


or 4r° 





1 
: 


and 
(i) Wa, T) = a’’9Q, 
(ii) Wr, 0) = 0, (1.2) 
(iii) lim W(r, T) = 0. 
Now we define the Laplace transform of W(r, 7) as 


W*(r, s) = [ Wer, Tew"? aT. 


) 


Thus Eq. (1.1) becomes 


27% 7% 
the solution of which is 

W*(r, s) = cl3,.(rs\’”) + DK3,2(rs'”’), 
where J;,. and K3,. are modified Bessel functions of the order 3/2. 

The boundary conditions now give 

c=0 and D = a’’Q/sK;,(as’”’). 
Therefore, 

W*(r, s) = a°’2K3,2(rs'””)/sK3,2(as””). 


From Laplace inversion formula 


3/2 c+i@ 1/2, 37 
7 = 3/2 oe a 2 K3,.(rs Je 
We, T) = Pelt, T) = “eK, (a8) ds. 





Jc-i@ 


For the order (n + 1/2), where n is an integer, the Bessel functions and modified 
Bessel functions both reduce to finite form. Thus we get 


'’ a’Q re**? (1 +78” ‘ 
wr, T) = Oars . ea, exp [(a — rs" Je . ds, 


2 [z, + (- - 1)r. | , 
rT a 


where J, and J, are the inverse Laplace transforms of 


s”' exp [(a — r)s'”] and s''?(a"' + 8'”)"' exp [(a — r)s'””] 


A reference to the Tables of integral transforms, edited by Erdélyi and associates [2] gives 


7 ,(r-—a _ r-—a r)| (‘= ) 
I= 1 ef (58) and I, = [1 — ext (Gat + : exp - ae ’ 
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where 


We thus have, 


71/2 _— 7 
vt, 1) = S24(t —1)] 1 — ent (Giant + 2) | ex ("S44 3) 
+1—- ef (Sa p72 a) (1.3) 


The value of w(r, T) is given by the expression (1.3) and satisfies the equations of 
motion as well as the boundary conditions. The transient part in (1.3), say K, can be 
written as 


2a*(r — a)Q (27)? ; e r ) 
ais °° 1 - Pla + a few ae 


T-_ a 
+ Pfr P| Ea 4} + oe 


1 sf : 
P(x) = an | exp (—4u°) du. 











where P(x) is the probability integral 


For large values of r, K « 1/r*°. Asr — ~, K approaches the value zero quite rapidly. 
K also approaches zero as T — ~, 
Now let a be the unit of length and 27’ = 1. Equation (1.4) then reduces to 


° 
K= “3 {2(r — 1)[1 — P(r)] exp ( — 4) + 1 — 2P(r — 1)}. 
The following table’ gives the values of K/Q for different values of r when 27’ = 1. 


Relation between r and K/Q. 








r K/2 r K/2 r K/2 r K/2 r K/2 
1.00 -—0.0000 1.60 -—0.0620 2.40 -—0.0495 3.60 -—0.0210 6.00 -—0.0046 
1.10 -—0.0227 1.70 —0.0628 2.60 -—0.0437 3.80  —0.0180 
1.20 —0.0338 1.80 —0.0624 2.80 -—0.0380 4.00 -—0.0154 
1.30 -—0.0485 1.909 -—0.0615 3.00 -—0.0324 4.50  -—0.0109 
1.40 -—0.0553 2.00 -—0.0598 3.20 -—0.0282 5.00 -—0.0079 
1.50 —0.0596 2.20 —0.0551 3.40 —0.0241 5.50 —0.0061 





The table shows that unsteadiness spreads from the surface of the sphere into the 
interior of the liquid and attains its maximum value for r = 1.70, approximately. There- 
after it gradually decreases and becomes inappreciable at a finite distance from the 





1The values of P(x) are taken from Biometrika tables for statisticians, vol. 1, ed. by E. S. Pearson 
and H. O. Hartley, 1956, pp. 104-106. 
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sphere. The region in the neighbourhood of r = 1.70 is consequently the region of maxi- 
mum unsteadiness. 
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DYNAMIC PROGRAMMING APPROACH TO OPTIMAL INVENTORY 
PROCESSES WITH DELAY IN DELIVERY* 


By RICHARD BELLMAN (The RAND Corporation, Santa Monica, California) 


Summary. The usual dynamic programming approach to inventory processes with 
delays in delivery leads to functions of many variables. This multi-dimensionality 
prevents the straightforward utilization of digital computers. 

Using a type of transformation previously applied in the study of engineering control 
processes, we show that a class of inventory processes with time lags can be treated in 
terms of sequences of functions of one variable, regardless of the length of the delay. 


1. Introduction. The problem of determining ordering policies which minimize 
the cost of operating supply depots and stockrooms is one which has attracted a great 
deal of attention in industrial and military circles in recent years. An analytic approach 
to these questions by way of functional equation techniques was inaugurated by Arrow, 
Harris, and Marschak, in a now classic paper, [1]. These investigations were extended by 
Dvoretzky, Kiefer, and Wolfowitz, [7], and Bellman, Glicksberg, and Gross, [6]; see 
also [2], and the books by Whitin, [8], and Arrow, Karlin and Scarf, [9]. 

Although this approach can be used to obtain analytic and computational solutions 
of a variety of processes in which there is no delay between an order for an additional 
supply of items and the delivery of these items, this method runs into dimensionality 
difficulties when time lags of more than a stage or two occur. If there is a delay of d 
stages in filling an order, the state of the system at any time is characterized not only 
by the present stock level, but also by the quantities on order which will arrive one, 
two, --- , d stages in the future. 

It thus appears that functions of d variables necessarily arise when point of regener- 
ation methods are employed to treat these processes. 

In several papers devoted to the study of control processes arising in the engineering 
world [3], [4], we have shown that in some fortunate situations certain preliminary trans- 
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formations enable us to bypass functions of large dimension. Problems which seemingly 
require functions of ten or twenty variables can actually be reduced to questions in- 
volving functions of one or two variables. This vast reduction in dimension renders a 
computational solution feasible. An example of this phenomenon in the field of scheduling 
is contained in our paper on the “caterer” problem, [5]. 

Here we shall consider the application of this type of transformation to inventory 
processes. The stochastic nature of the process introduces certain features of compli- 
cation. Although under appropriate assumptions of linearity exact solutions can be 
obtained, we shall concentrate here only on making the problem susceptible to compu- 


tational solution. 


2. Description of process. In this simplified model of an inventory process, we 
consider the question of stocking a single item which is subject to a stochastic demand 
at certain specified times which we denote by 0, 1, --- , NV. Since there is a penalty at- 
tached to inability to supply these demands, quantities of this item are purchased 
throughout the process at these same time points. This outage or penalty cost is a 
prescribed function of the excess of demand over supply. We suppose that there is a time 
lag of d units between an item being ordered and its delivery. 

In many cases, there is an additional cost incurred in storing items. Since this 
introduces no additional analytic difficulties, but does complicate the arithmetic, we 
shall ignore this cost in our future discussion. 

Assuming that we are given the purchase cost, the outage cost, and the distribution 
of demand at any stage, we wish to determine ordering policies which minimize the 
total expected cost. 

Complicated as this problem is, it is one level simpler than many which arise in 
practice in which the distribution of demand is not completely known. These will be 
discussed at another time; see Dvoretzky, Kiefer and Wolfowitz [7]. 


3. Analytic formulation. A discussion of a straightforward dynamic programming 
formulation of processes of the type described in the foregoing paragraph may be found 
in [2]. As mentioned above, this direct approach leads to dimensionality difficulties. 
Here we postpone the introduction of dynamic programming techniques until an appropri- 


ate point. 
Let us define the following quantities: 
a. x, = the stock level at the nth stage, prior to the delivery of the quantity ordered 


at the nth stage and the demand at the nth stage. 
b. y, = the quantity ordered at the nth stage, 


n=0,1,---,N —d; (1) 
c. z, = the demand at the nth stage, n = 0,1, --- , N; 
d. w, = the actual cost of the first n stages, n = 0, 1, --- , N. 
All four of these quantities are stochastic with the exception of x» and y». The quantity 
y, depends upon the observed stock level and the quantities ordered at preceding stages, 
Yn-1 5 °** » Yn-a - It is assumed here that y, must be determined before z, is observed. It 
will be clear that either case could be discussed. 
The following relations are then immediate, 


Lart = In — Zn + Yara » (2) 
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where y, = Ofork < Oandk > N — d, and 
W413 = Ww, + G(Yn+1) + O(2n41), y= 0, :. htt (3) 


Here g(y) measures the cost of ordering y items and ¢(x) is the cost of an outage of x 
items. A negative stock level is equivalent to outage. 
It follows that 


Wy = bs g(yi) + > g(x), (4) 


where the x; are given by 
X= MX — (tat: + 2-1), a2=1,2,--- ,d— Il, (5) 
Li, = Xo — Bo Hi + ++> $2-1+) + Pot Ht -°> + H-<; 
a=d,--- ,N. 
Substituting, we see that the total cost of an N-stage process may be written 


N-—d 


ws = 2 ays) + ols — (@o +a + +++ + %-1)) 
_ _ + (to — (o +21 + o°** + 20-1) + Yo) 
+ (to — (oo 2+ e+: + z.) + Yo + %) 
+ +++ + Olt — (% fai $+ °°* + en-1) 
+nth > +--+ + w-2- 


We now wish to determine the policy functions, {y;}, stochastic functions of the 
z; , which will minimize the expected total cost. This expected value is taken over the 
random variables, {z,;}, independent random variables with a known distribution 


(6) 


function. 


4. Dynamic programming formulation. It is clear that in seeking the policy functions 
which minimize Exp (wy), we can ignore the term 


Exp | D> d(xi — (eo tees + 2) | 


which is independent of the y; . 
Let us then, for V > d, introduce the sequence of functions defined by 


N-d 


N-1 
fy(a) = Min Exp b> g(y;) + + b(% — (2 +2; +--+ +e) 
P k=0 ked—-1 


(1) 
+ytwte:: + new) } 


The notation Min, indicates that the minimization is over all policy functions. 
We have 
f(z) = Min Exp [g(yo) + O(t%o — (@o +21 + +++ + 24-1) + Yo). (2) 


vo Zo.ta. 1 


To obtain a recurrence relation connecting fy(2) and fy_,(x), observe that the effect 











402 RICHARD BELLMAN (Vol. XVIII, No. 4 


of an initial choice of 7» and an initial demand z is to produce a similar sum with 2» 
replaced by % + Yo — 2 . Taking account of the initial terms g(yo) + Exp,, ¢(%> + 
Yo — Zo), we expect to obtain the following recurrence relation upon using the principle 


of optimality, [2]. 


f(x) = Min Ee + Exp G(X — (2p + ++ + 24-1) + Yo) 
ve (3) 


7) 


fu-1(Xo Scio Zo + y) aH, | ? 


70 
N=d+1,::: 
5. Rigorous derivation. In order to derive this recurrence relation in a rigorous 
fashion, we must make precise what is meant by a sequence of ordering functions and 
what is meant by minimization over these functions. 


The quantity y; is a stochastic variable, dependent upon 2 , 2% , 2%, °** , 2:-, and 
Yo, 41, °°* »Yi-1,, TY: = Yr (Zo, 2 9° °° » Benny Yo, Yt» ° °° » Yo- ; x). The multistage 
nature of the process shows that fy(*) may be written 

fy(z) = Min [Exp[ Min  [Exp[---]]]]. (1) 
vo (zo) zo ¥i(z0,¥0;Fo) 21 


It follows that we can write 


fv(z) = Min [g(yo) + Exp (2% — (o + 21 + *** 4+ 2-1) + Yo) 


eat (2) 
+Exp Min _ [Exp[---]]]. 
Zo Va (20,0320) Z1 
The expression 
gv-1 = Exp Min _ ([Exp[---]] (3) 
Zo «= (20,0320) z1 


looks like Exp,, fy-1 (vo + Yo — 20), with the difference that the minimization in (1) is 
over the functions ¥o(%o), ¥: (20 , Yo ; Xo), and so on, while the minimization in (3) is over 
the functions y,; (20 , Yo } Zo), Y2 (20 , 21 » Yo» Yi j Lo), and so on. Observe, however, in 
what fashion the variables zp , yp and 2» enter in (3). They occur only in the combination 
Zo + Yo — 2% . Consequently, minimization over ¥; (Zo , Yo ; Zo) is equivalent to minimi- 
zation over ¥; (%o + Yo — 20); similarly, minimization over yz (20 , 21 , Yo» Yi j Lo) 18 
equivalent to minimization over Y2 (2; , ¥1 } 2o + Yo — 20), and so on. 
It follows that 
gv-1 = Exp fy-s(%0 + Yo — 2)- (4) 
Ze 


This establishes the formula of (4.3). 


6. Discussion. We see then that the solution of the original inventory problem 
can be treated in terms of a sequence of functions of one variable, regardless of the size 
of d. Under various simplifying assumptions concerning purchasing cost and outage 
cost, we can determine the structure of the optimal policy; see, for example, the dis- 
cussion in [2], [9]. 

If, however, we take into account fixed costs which occur in realistic processes, 
formidable analytic difficulties arise. Nonetheless, the computational solution can 
readily be obtained using (4.3). 
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AN ELEMENTARY DISCUSSION OF DEFINITIONS 
OF STRESS RATE* 


By WILLIAM PRAGER (Brown University) 


1. Introduction. The simplest constitutive equation considered in the theory of 
plasticity only involves the tensors of stress and rate of deformation and described 
rigid, perfectly plastic behavior. When elastic effects are to be included in the analysis, 
this equation is assumed to apply to the plastic part of the rate of deformation, to which 
an elastic part must be added before the total rate of deformation is obtained. In a 
similar manner, a simple constitutive equation for a viscoelastic material can be estab- 
lished by adding an elastic rate of deformation to the rate of deformation of a viscous 
fluid. In both cases the elastic rate of deformation is usually written as a function of an 
appropriately defined rate of stress. 

This stress rate must obviously satisfy the following condition: if a stressed con- 
tinuum performs a rigid body motion and the stress field is independent of time when 
referred to a coordinate system that participates in this motion, the stress rate vanishes 
identically. As is readily seen, this restriction is not severe enough to lead to a unique 
definition of stress rate. Indeed, from one definition that satisfies this condition another 
one may be obtained by adding terms that contain the rate of deformation. In a rigid 
body motion, these terms vanish, and the second definition reduces to the first. Since 
the condition imposed on the definition of stress rate only concerns rigid body motions, 
it will be satisfied by the second definition if it is satisfied by the first. As a consequence 
of this freedom of choice, many definitions of stress rate are found in the literature. 

A similarly embarrassing choice offers itself when one attempts to define finite strain. 
Potentially, any tensor formed from the displacement gradients qualifies as strain 
tensor if it vanishes identically for all rigid body motions. Depending on the field of 


*Received Feb. 6, 1960. This paper is based on a report distributed under Contract Nonr 562(10) 
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application, one or the other definition of strain may be more convenient to use. Whereas 
the relative merits of the various definitions of strain are well understood, those of the 
various definitions of stress rate have not been discussed in similar detail. Moreover, 
some of these definitions are usually presented in a manner that presupposes familiarity 
with general tensor calculus. In this paper, only Cartesian tensors are used; the principal 
definitions of stress rate found in the literature are derived in an intuitive manner, and 
their suitability for the mathematical description of elastic, plastic behavior is discussed. 

2. Notations and basic relations. Let ¢ denote the time and 2,(i = 1, 2, 3) a fixed 
system of rectangular Cartesian coordinates, and denote partial differentiation with 
respect to these independent variables by the operators 0, and 0; , respectively. 

A tensor of the second order, 7';; , can be written as the sum of its symmetric and 


antisymmetric parts, which will be denoted as follows: 


Tis) = 27; + T;,), 
(1) 
T ii = a(7';; ans T ji). 
If v;(%; , X2 , Xs , t) is the time-dependent velocity field of a continuum, the rate of 


deformation is defined as the symmetric part 0,,v;, of the tensor 0,v; , and the rate of 
rotation as its antisymmetric part 0;,v;; 

The acceleration of a typical particle is the material derivative v; of its velocity, 
which is defined by 


Vv; = Ov; + 0; 0); , (2) 


where the repeated subscript indicates summation in the usual manner. The material 
derivative of the infinitesimal vector ds; joining adjacent particles is 
7) és / 
(ds;) = 0,0; ds; . (3) 
Because it will suggest a possible definition of stress rate by analogy, we consider 
the second rate of deformation. If ds; and 6s; are the infinitesimal vectors joining the 
typical particle to two neighboring particles, it follows readily from (3) that the material 
derivative of their scalar product is 
(ds; 68;) = 2 0:,%;) ds; 68; . (4) 
Forming the material derivative of (4), we obtain, after some simplifications, 
(ds, 6s;)° = 2[(0..0;)) + Ov, Auv;) + O% Ou0;)] ds; 6s; . (5) 
Since the coefficient of 2ds;6s; in (4) is called the rate of deformation, the corresponding 


coefficient in (5), i.e. the contents of the bracket, may well be called the second rate of 
deformation. Higher rates of deformation have been discussed by Rivlin and Ericksen 


(1]*. 


The oriented surface element formed by the infinitesimal vectors ds; and 6s; is defined as 
dA; oS Gy ds; 6S; ; (6) 
where ¢; ;, is the alternating tensor. It is readily verified by means of (3), that the material 


*Numbers in brackets refer to the Bibliography at the end of the paper. 
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derivative of this surface element can be written as 
(dA,)° = 0; dA; — 0; dA; (7) 


(see, for instance, [2, p. 598]). 
3. Jaumann’s definition. If o;;(x,, x2, x3, t) denotes the time-dependent stress 
field in the continuum, the material derivative 


0; ; = O09 5; + Vv; Ono; ; (8) 


indicates the time rate of change of a typical stress component at a particle of the con- 
tinuum, the stress tensor being referred to the fixed coordinate system. 

When the continuum and its stress field perform a rigid body motion other than a 
mere translation, ¢,; does not vanish identically. Accordingly, (8) is not acceptable as 
a definition of the stress rate for use in constitutive equations. A suitable definition is 
obtained as follows. 

Let P, Q, , Q2 , Qs be neighboring particles such that the lines PQ, (a = 1, 2, 3) 
indicate a system of principal axes of the rate of deformation at the particle P and the 
time ¢, and denote the positions of these particles at the time t + dt by P’, Q! , Q3 , Q3 .. 
As the shear rates for the principal axes of the rate of deformation vanish, the lines 
P’Q! are orthogonal. The stress rate could be so defined that it vanishes if the stress 
tensor at P has the same components with respect to the axes PQ, at the time ¢ as the 
stress tensor at P’ has with respect to the axes P’Q/ at the time ¢ + dl. 

Let A* be the unit vector of the direction PQ, and A$ + (A%)°dt that of the direction 
P’Q, . Since the axes PQ, have the rate of rotation @ ;,v;,; , we have 


(AS) = OF - (9) 


Now, a typical component of the stress tensor o,;; with respect to the axes PQ, is given 
by the expression o;;A°A* . According to (9), the material derivative of this expression is 


(ao; r*¥)° = (a;; -+- Cx; 0 iV] + Ok 0 ey ° (10) 


If the material derivatives of all these stress components are to vanish, the parenthesis 
in (10) must vanish. The expression in the parenthesis thus constitutes a suitable defi- 
nition of stress rate. On account of the symmetry of the stress tensor and the anti- 
symmetry of the rate of rotation tensor, this definition, which is due to Jaumann [3], 
may be written as 

of, = 04; — On OV;) — Tie Away - (11) 
Jaumann’s work does not seem to be well known: the definition (11) is frequently used 
in the recent literature [4] without reference to Jaumann. 

Since Jaumann’s stress rate o/; measures the rate of change of the stress components 
with respect to a rectangular Cartesian system that participates in the rotation of the 
material, ¢/; = 0 implies that the invariants of the stress tensor are stationary. 

4. The definition of Cotter and Rivlin. Let /* (2 = 1, 2, 3) be three linearly inde- 
pendent vectors emanating from the typical particle P at the time ¢. During the time 
interval (¢, ¢ + dt), the neighborhood of P undergoes the infinitesimal affine transfor- 
mation specified by tensor 0,v, . If the vectors /¢ are subjected to the same transformation, 


ve have, in analogy to (3), 


(12)” = 3,12. (12) 
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When the continuum and its stress field perform a rigid body motion, the material 
derivative 
(0, ;155)° = (65; + on; OM + ou dv, )LiG (13) 
vanishes. In view of the symmetry of the stress tensor, this remark yields the following 
definition of stress rate, which is due to Cotter and Rivlin [5]: 
oi; = Oi; + OMeon; + INK: - (14) 
Note that the right-hand side of (14) has the same structure as the expression in the 
bracket in (5) which can be written as 
(QWs) = (8¢0)° + Owe Oud;) + OM, Aun - (15) 


Cotter and Rivlin [5] have shown how higher rates of stress can be obtained by a 
continuation of the procedure that led to (14). While the definition (14) lends itself 
more readily to this generalization than the definition (11), it has the drawback that 
a’! = 0 does not imply stationary behavior of the stress invariants. 

5. Oldroyd’s definition. Let the stress tensor at the particle P and the time ¢ be 
written as a linear combination of the dyadic products of the vectors 7 : 


Ci; = CaplSlh . (16) 
On account of (12), the coefficients c,, of this combination will be stationary if 
01; = Sag OU 1G + 8a 0, 1ilt (17) 
On; OV; + O%~ OY; 


In view of the symmetry of the stress tensor, it follows from this remark that 


€ a . 
of! = O45 — O% OY; — Ojn OW; (18) 


is an acceptable definition of stress rate. 


This definition, which is due to Oldroyd [6], has the same disadvantage as that of 
Cotter and Rivlin: the stress invariants need not be stationary if o//’ vanishes. 


6. Truesdell’s definition. The infinitesimal force transmitted across the oriented 
surface element dA; is 
dP; = a;; dA; . (19) 
Let this force be written as a linear combination of the vectors /< : 
dP; = dC, I; . (20) 
According to (12), the coefficients dC, of this combination will be stationary if 
(dP;)° = dC, dv;le = dP, 0; . (21) 
Substitution of (20) into (21) and use of (7) furnishes o//’’"dA; = 0, where 
OL" = a3; + 04; Ove — Tin OW; — oj, OY; - (22) 


As follows from its derivation, Eq. (22) represents an acceptable definition of stress rate. 

This definition, which is due to Truesdell [2, p. 604], may seem more far-fetched 
than the preceding ones, because it is based on the consideration of elementary forces 
rather than stresses. It can be shown, however, that Truesdell’s stress rate is closely 
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related to the partial derivative of Kirchhoff’s stress tensor with respect to time, the 
other independent variables being the rectangular Cartesian coordinates of the particles 
in some reference state (see [7], Chap. IX, Sec. 4). Note that the vanishing of Truesdell’s 
stress rate does not imply stationary behavior of the stress invariants. 

7. Use in constitutive equations of plasticity. In the theory of perfectly plastic 
solids, a state of stress is said to be at or below the yield limit according to whether a 
certain function of the stress invariants, the yield function, is zero or negative. States 
of stress that furnish positive values of the yield function cannot be supported by these 
solids. 

The constitutive equations describing elastic, perfectly plastic behavior are generally 
conceived as resulting from the superposition of the rates of deformation of an elastic 
and a rigid, perfectly plastic constituent. To facilitate this superposition, the consti- 
tutive equation of the first constituent is usually assumed to establish a one-to-one 
correspondence between the rate of deformation and the stress rate. The second con- 
stituent is supposed to be rigid under stresses below the yield limit, and to flow plastically 
under stresses at the yield limit. 

In the composite elastic, perfectly plastic solid two criteria are therefore used to 
judge a given variation of stress in time: the stress rate in the elastic constituent, and 
the rate of change of the yield function in the plastic constituent. To avoid contradictions, 
the yield function should be stationary when the stress rate vanishes. A similar condition 
arises in the theory of elastic, work-hardening solids, where a vanishing stress rate 
should go with a stationary state of hardening. Only Jaumann’s definition of stress 
rate satisfies these conditions. For use in the constitutive equations of plasticity, Jau- 
mann’s definition is therefore preferable to the other definitions of stress rate discussed 
in this paper. 
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ON THE FLOW OF A VISCOUS ELECTRICALLY CONDUCTING FLUID* 
BY H. P. GREENSPAN (Massachusetts Institute of Technology) 


The equations governing the steady two-dimensional flow of a conducting fluid 
whose free stream direction is parallel to an applied magnetic field can be linearized, 
by an extension of the Oseen treatment of viscous flows, and reduced to the dimension- 


less form 


dy ,. 24) 
a(ay ae * Fag] = % () 
4— (9A _ av) - . 
AA ( a = ©. (2) 


Here the stream function, ¥(z, y) and vector magnetic potential A(x, y) are defined by 
the relations 
q= V X Wi, y)is 
H= V X A(z, y)is - 
In addition, 8 = (uH>)/pV>), and « = ov. (H, and V, are the free stream values of 
magnetic intensity and fluid velocity; the other symbols have their usual meanings.). 
The details of the reduction can be found in Ref. [1]. 
A transform analysis of the flow past a rigid impermeable semi-infinite flat plate, 
[1], results in a solution of the foregoing equations which is of the form 
¥ 
A 


II 


Ef(n), (3) 


Eg(n), 


with 
(— + in)? = xt ty. 


Equations (1) and (2) admit general solutions of the form expressed in (3) and an 
entire class of parabolic flow problems is then explicitly soluble in terms of known 
elementary functions. Typical among these would be the flow in an aligned magnetic 
field past a parabolic cylinder which may be injecting another conducting fluid into 
the main stream. At least three different regions and possibly five independent param- 
eters are involved in this complex interaction and an exact numerical solution of the 
non-linear equations would require an immense amount of work. The linear theory, 
however, affords an easy way of examining the qualitative effects of a variation of any 
particular parameter on the fluid motion or magnetic field because explicit analytical 
solutions are readily obtainable. In other words, a rapid qualitative survey of any 
particular aspect of the phenomenon can easily be made using the linear theory and 
this can then be followed by a quantitative analysis if the results warrant a more exact 


investigation. 





*Received March 23, 1960. This research was partially supported by the U. S. Air Force under 
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‘ 
4 


The substitution of the functional forms of Eqs. (3) into (1) and (2) yield 
—4 1 | f”’ ey 

773 ae, -S [5 + apt + 09 - a | 

G+) dnl 2 


t a’ ”, (4) 
+ptaG[S-14+0 + a9 — | = 0 


and 
3g’ + ef — nf’ — 9g + ng’) = 0. (5) 


The former equation is satisfied (fortuitously) if 


l :' 
ay Ul’ — f+ af’ + Bg — 19')] = 0 


af’’ — } + af’ + Bg — ng’) = K, (6) 
where K is a constant. The general solution of Eqs. (5) and (6) is, for 8 < 1 (super- 


Alfvén flow) 


‘ G(Qn) K yy 

f(n) + ag(n) = aon + a, - AS) —_ a3 , (7) 
G ; 

f(n) + aog(n) = bon + by ee... . , (8) 
w w 


where 


G(n) = 7 erf » + 


? 


exp (—17') 
r' /2 


OV =1l+e—d, W=-ltet+aA, A=(1—&)? + 48)", Za» =1l—etA, 


and d> , a, , 0) , b, , K are all arbitrary constants. If there is more than one parabolic 
domain i.e. body plus fluid, the general solutions appropriate in each region must be 
joined by matching conditions at the boundaries. 

As a particular example consider the flow past a semi-infinite flat plate along which 
the horizontal velocity component is zero and the vertical component is —d ¥(x, 0)/dx = 
— V/2z*, x > 0 (the constant V is negative for injection and positive for suction). The 
boundary conditions require that {(0) = V,g(0) = f’(0) = Oand f’(@) = 2. The solution 


is then 





so r “9 — F 2 = Q 
f(n) = 2n + | 20a + 77 VO ~-e ea |, — ans) He) 


ae ‘ ©) 
_ Vig? — «? 
-= | 200, + nl? i = asa, [ou - aw)" Hen ’ 
1 2 
g(n) = 29 — (Qn + 2? V(2? — w")ay)(a,2 — ago)” Aon 
(10) 


+ (22 + 2'?V(2? — w’)a2)(a,2 — aw)" 


’ 


H(wn) 
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Fig. 1. The flow past a flat plate which is injecting fluid into the mainstream. Solid curves represent 
streamlines and dashed curves are field lines. The injected fluid is the same as that of the free stream. 


where 
H(n) = nerfen +a (1 — exp — 7’). 

A typical injection flow pattern is shown in Fig. 1. The two distinct flow regimes which 
occur are separated by the zero streamline y = éf(n) = 0. (That is,& = Oandz+ 4° = 
1y°4-° where f(7) = 0.) The interior domain is occupied by the injected fluid and the 
exterior realm consists of the fluid from the free stream. Two distinct magnetic fields 
also occur; an induced field occupies the domain interior to ég(n) = O and the applied 
field occupies the exterior region. Descriptively speaking, the free stream lines are pushed 
outward by the emerging fluid and drag the magnetic field lines behind them. The 
induced interior magnetic field closely resembles the field of a semi-infinite solenoid of 
current. 

The solution of the more general problem involving the flow past a parabolic cylinder 
(which may be injecting or drawing fluid or impermeable) can be solved in as straight- 
forward a manner as the preceding example. The general solutions, valid in each region, 
must be properly joined at the interface boundaries. Across the parabolic interface 
separating two fluids it is required that the fluid velocity, viscous stress, tangential 
component of magnetic intensity, V x Ai, , and normal component of magnetic field 
nV X Ai; be continuous. (The latter two conditions hold across any surface.) This 
corresponds to requiring the continuity of f, f’, f’, ug, and g’. The tangential velocity 
component at the surface of a body is zero although the normal component may be 
prescribed to simulate injection or suction. (It is necessary, however, to maintain a 
parabolic flow.) The effects of magnetization can also be examined by prescribing a 
non-zero boundary condition on g(0). The arbitrary constants are easily evaluated and 
the solution is explicitly determined in terms of elementary functions, the most compli- 
cated of which is erf n. The location of an interface such as that between injected fluid 
and free stream is obtained from the condition that it constitutes part of the zero stream 
line éf(n) = 0. 


A particular solution, worth quoting because of its simplicity, is that for the perfectly 
x 


conducting flow past the impermeable parabolic cylinder 7 = 7 


f(n) = g(n) = 2y — 9") - 20 = 8) ta — 9") — HA — 8)""9")] 
; erfe (1 — B)'’*n* ; 





The first term of the right hand side represents the inviscid solution and the second is 
the viscous magnetohydrodynamic correction. 
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The linear theory is, of course, grossly inaccurate near stagnation points of either 
fluid or magnetic field. Other sources of extreme errors may arise because the actual 
magnetic force is replaced by its component normal to the free stream direction. This, 
for infinite conductivities, requires shear discontinuities to maintain the dynamic 
equilibrium of current carrying interfaces. It should also be noted that the linearized 
version of the infinite conductivity condition q X H = 0 need not require the alignment 
of field and fluid if they are not already so directed somewhere within the perfectly 
conducting region. If the linear theory is not applied to problems for which it is obviously 
inappropriate, fairly good qualitative results can be expected. 

Substantially correct quantative results may also be obtainable by rescaling the for- 
mulae developed from this theory in accordance with the modification introduced by 


Lewis and Carrier [2]. 
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BOOK REVIEWS 


(Continued from p. 394) 


use illustrated in the calculation of definite integrals. The third chapter, entitled ‘““Les Formules Ré- 
ciproques de Hankel,’”’ is misnamed since it contains not, as we should suppose, a derivation of the 
Hankel inversion theorem, but an account of the formal properties of Hankel transforms (the main 
theorem being stated without proof). The next chapter is concerned with applications of these trans- 
forms to the solution of certain simple boundary value problems—Laplace’s equation for a wedge by the 
Mellin transform, Laplace’s equation and the diffusion equation by the Hankel transform. All of these 
problems have been discussed in earlier books by Titchmarsh, Tranter and the present reviewer. The 
book ends with a brief discussion of the solution of the classic pair of ‘“Titchmarsh’’ dual integral 
equations. ‘ 

The book is written clearly and provides for the French-speaking world an admirable introduction 
to these transforms. Anyone whose native language is English would, however, be better advised to 
read C. J. Tranter’s “Integral Transforms in Mathematical Physics’ (Wiley, New York, 1951). 


I. N. SNEDDON 


Incompressible aerodynamics. Editor Bryan Thwaites. Oxford University Press, New 
York, 1960. xx + 636 pp. $12.00. 


The present book is an account of the theory and observation of the steady flow of incompressible 
fluid past airfoils, wings, and other bodies of aeronautical interest. Written by sixteen collaborators, 
it is the first of a new series sponsored by the Aeronautical Research Council to succeed the well-known 
pair of works edited by Sydney Goldstein and by Leslie Howarth. The present volume, together with 
a second volume on laminar boundary layers and a third on turbulence yet to appear, will recover, in 
the main, the topics in the Goldstein volumes. 

Following an inauspicious introductory chapter, the second chapter presents a brief modern dis- 
cussion, from a utilitarian point of view, of the computation of chiefly turbulent boundary-layer charac- 
teristics, including the results of D. Coles and F. Clauser. The next chapter compares the results given 
by theoretical models of real flows with those observed. Delightful photographs here illustrate the need 
for fast exposures to depict what actually happens, for example, in the wake of a disk, the more familiar 
pictures being revealed as smears of the phenomena. The usual two-dimensional airfoil theories are 
briefly but nicely set out in Chapter IV. Lighthill’s theory for rendering perturbation results uniformly 
valid at stagnation points is included, as is necessary in a modern treatment, and the chapter goes on 
to treat other recent results. The pristine theories of separated flow remind us (and perhaps we need it) 
that the two-dimensional, incompressible flow past airfoils is still a living, vital subject, providing, in 
fact, a framework suitable for the study of separated flow. The next chapter on the uniform viscous flow 
past airfoils includes a critical generalized view of the Kutta condition; it deals with its subject from 
the point of view (maintained throughout the volume) of successive approximation to the actual vis- 
cous flow via the first inviscid approximation, the first viscous approximation (the boundary layer), the 
second inviscid approximation (displacement effects), the second viscous approximation, etc. 

Following two interesting chapters on boundary-layer control and the thickness effects of three- 
dimensional wings, the subject of three-dimensional lifting wings, probably the central chapter in such a 
book, is treated. This is introduced with some remarkable idealized sketches of the separated flow pat- 
tern and associated trailing vortex system produced by a flat, circular disk at varying incidences. These 
lead to the Lanchester-Prandtl vortex systems and to the recent Roy model for the separated delta 
wing. The classical linearized lifting-surface theory and slender-body theory are then presented, empha- 
sizing more recent developments; this material overlaps somewhat with Volume VII of the Princeton 
Series. Chapter IX on bodies of revolution follows the same pattern. Wing-wing and wing-body inter- 
ference is treated in the next Chapter. The discussion of the importance of nonlinear effects (in the 
boundary conditions) near the junction of two thin rectangular wings in cruciform arrangement is very 
interesting. Incompressible flow in axial-flow compressors, a pleasant surprise in a book of this type, 
is next considered, most of the chapter being devoted to actuator-disk theory. Concluding Chapter XII 
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is concerned with miscellaneous fringe topics, notably recent and important jet flap research and shear 
flows past simple bodies. 

This thoroughly modern treatment is distinguished for its attention to observation, its revealing 
three-dimensional sketches and photographs, its carefully reasoned formulation of idealized models, 
and its ever present experimental verification. If, indeed, the book has a categoric shortcoming, it is 
in an incomplete interpretation of the results of the analysis itself in respect to the feedback on the 
a priori physical reasoning and to the deeper connections between related analyses. As an example, 
consider the leading-edge singularities given by linearized airfoil theory. In accord with the point of 
view of the book, the formal manner in which the Lighthill theory for uniform validity manages to tuck 
the leading-edge singularity just the proper distance inside the wing by expanding the independent 
variables deserves interpretative discussion. For various cases, Goldstein and Van Dyke (dubbed Dyke 
in the book) have proposed multiplying the linearized solution by a factor which is the ratio of the 
deduced correct behavior at the leading edge to the linear behavior. This does not change the behavior 
of the solution elsewhere, within the approximation, and leads, as is mentioned in the book, to the 
same correct result as the Lighthill theory at the singularity. One wishes to see established the con- 
nection between the two diverse approaches, especially considering the success Lighthill has had in 
providing similar physical interpretations of his theory for cases of supersonic wave propagation (volume 
VI of the Princeton series). 

Didactically viewed, it isn’t likely that the book will be used as a text, but it will serve well to flesh 
out with supplementary reading briefer treatments of the subject. The work can be read by a first- 
year graduate student without difficulty or need of classroom support. At the other extreme from the 
Princeton Series, the chapters are so thoroughly integrated and otherwise worked over by the editor 
that, like the Goldstein volumes, the names of the contributors are not associated with their original 
contributions. Attributed articles carry the concomitant features of individual recognition and re- 
sponsibility, whereas the need to coordinate the efforts of collaborators is also obvious. One wonders 
whether the integrated but signed chapters of the Howarth volumes do not represent the best compro- 


mise solution to the problem of producing compendia. 
JosePH H., CLARKE 


Statistics. By Earl K. Bowen. Richard D. Irwin, Inc., Homewood, Illinois, 1960. xi + 
415 pp. $6.50. 


This text represents an introductory course for students of business administration, industrial 
management, and economics. Realistic problems are integrated with a systematic presentation of 
principles and methods and no formal training in higher mathematics is needed for understanding it. 


There are many questions and problems illustrating the material of each chapter and an extensive bibli- 


ography. The book should make interesting and useful reading for anyone wishing to acquaint himself 
with the multitudinous uses to which statistics can be put in the context of economic problems. 


WALTER F. FREIBERGER 


Mathematical methods of operations research. By Thomas L. Saaty. McGraw-Hill Book 
Co., New York, Toronto, London, 1959. xi + 421 pp. $10.00. 


After a preliminary chapter entitled ‘‘History and Concepts,’’ the book is divided into four parts. 
The first of these contains a chapter discussing scientific method in operations research, a chapter on 
mathematical existence proofs and symbolic logic, and a chapter summarizing, with some operations 
research examples, the use in problem formulation and solution of such mathematical techniques as 
differential equations (Lanchester’s application to combat), generating functions (probability), integral 
equations (renewal problems), and others. 

Part Two, entitled “Optimization, Programming and Game Theory,’ contains a chapter each on 
these subjects. The first of these is a collection (65 pages) of mathematical methods which the operations 
researcher should welcome, including the solution of inequalities, the use of Langrange multipliers, and 
the calculus of variations with inequality constraints. A chapter on linear and quadratic programming 
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is followed by a short chapter on the theory of games. A special section is appended to Part Two listing 
some mathematical exercises. 

Part Three contains a chapter on probability and another treating some applications in inventory 
processes and system reliability. There follows a review chapter on the elements of classical statistics, 
and finally the resume of queueing theory which the author published in Operations Research (vol. 5, 
no. 2, April 1957, pp. 161-200). A set of exercises is appended. 

In Part Four, an essay chapter entitled “Some Thoughts on Creativity”’ is followed by a collection 
of 95 problems, These problems include mathematical and logical teasers of the sort popular in scientific 
journals in recent years, as well as some suggested projects in operations research. 

The chapter-by-chapter bibliographies are well done. The interest of the books is, as the title im- 
plies, in mathematical methods, and it should be found useful by the applied mathematician seeking a 


collection of this size of topics relevant to operations research. 
HERBERT P. GALLIHER 


Turbulent flows and heat transfer. Editor: C. C. Lin. Princeton University Press, New 
Jersey, 1959. xv + 549 pp. $15.00. 


This book is intended to deal with the problems of turbulent flow and conductive, convective, and 
radiative heat transfer. Certainly many of the articles appearing in the volume attest to the high 
standard of excellence which may be expected from some of the outstanding contributors to it. On 
this basis alone it is to be recommended to those people who wish to obtain an appreciation of the sub- 
jects treated, particularly of the complex area of turbulence and transition to turbulence. 

As a connected work the volume falls far short of what might have been expected. This is particu- 
larly true of those sections concerned with heat transfer. The separate sections in the book span a range 
which varies from a handbook engineering style to highly sophisticated theoretical discussions. The 
caliber of the articles has an equally wide spread going from poor to excellent. Furthermore, as a result 
of the disconnected presentation there results considerable duplication from different points of view 
and as a consequence the reader is oftimes left with a confused and, in some cases, erroneous impression 
of the field. 

The nine articles in the book in their order of presentation are: 

A. Transition from laminar to turbulent flow, Hugh L. Dryden; B. Turbulent flow, Galen B. 
Schubauer and C. M. Tchen; C. Statistical theories of turbulence, C. C. Lin; D. Conduction of heat, 
M. Yachter and E. Mayer; E. Convective heat transfer and friction in flow of liquids, R. G. Deissler 
and R. H. Sabersky; F. Convective heat transfer in gases, E. R. van Driest; G. Cooling by protective 
fluid films, S. W. Yuan; H. Physical basis of thermal radiation, S. S. Penner; I. Engineering calculations 
of radiant heat exchange, Hoyt C. Hottel. 

It is in Sections A to C that this reviewer found a sensible and orderly presentation of the material. 
Section A opens the volume with a discussion of transition from laminar to turbulent flows essentially 
from an experimental point of view, and is followed in B with a discussion of turbulent shear flow from 
a semi-empirical and experimental point of view, while in C turbulence is discussed from the purely 
deductive statistical point of view. 

Dryden’s contribution in A is an excellent compilation of all of the presently known experimental 
facts on the subject of transition to turbulence, and is a fine presentation of how various flow factors 
influence transition in both incompressible and compressible flows. This reviewer would have liked to 
have seen more detailed considerations on the three-dimensional nature of the transition process and 
perhaps brief outlines of some of the major existing theories on stability and transition to turbulence. 
However, this in no way detracts from the value of this section. 

The article by Schubauer and Tchen in B emphasizes the mean properties of turbulent shear flows 
without discussing the structure of turbulence. Here is to be found an excellent description of the 
various phenomenological theories of turbulent flow and turbulent boundary layers. A good discussion 
is given of the wall law and velocity defect law, Clauser’s equilibrium boundary layer concept, and 
Coles’ law of the wake. This reviewer was pleased to see the authors’ statement in connection with so- 
called turbulent boundary layer skin friction “theories’’ that ‘‘some empirical form of skin friction 
coefficient may often be more useful in practice.” 

To the reader unacquainted with statistical theories of turbulence, C, C. Lin’s article in C will indeed 
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be a welcome addition. The section restricts itself primarily to the case of homogeneous turbulence and 
the subject is presented from the author’s own unique point of view. A clear discussion may be found 
here of velocity correlations and the spectral theories of turbulence. These methods are then applied 
to the study of the nature of turbulent motions and similarity considerations are presented. Turbulent 
motion in a compressible fluid is briefly discussed. 

The remaining sections of the book deal with conductive, convective, and radiative heat transfer 
with special emphasis on high speed flows. The longest section of this part which is intended to analyze 
convective heat transfer is Section F by van Driest. Unfortunately many of the topics contained therein 
are unrelated to the topic at hand, so that some of what was presented in A and B on transition and 
turbulent flows is repeated here in a much less satisfactory manner. The remainder of the section appears 
to deal almost exclusively with particular problems the author himself has been connected with. After 
reading this section this reviewer was not left with any clear picture or understanding of the general 
area of convective heat transfer in high speed flows. The shortest section in the volume is H by S. S. 
Penner which trys to present in 12 pages the basic ideas of thermal radiation. In this task the author 
is not successful and the uninitiated reader will find it difficult to actually get a clear picture of the prob- 
lems involved, particularly in relation to the calculation of gas emissivities. 

Section D by Yachter and Mayer gives a brief discussion of conduction in solids covering a variety 
of topics from a formal mathematical discussion of linear heat conduction problems to practical con- 
siderations in minimum weight design. Isolated problem solutions are given and brief remarks of interest 
to those familiar with the problems involved are given. Deissler in E discusses turbulent heat transfer 
in liquids from a semi-empirical and experimental point of view with no reference to the discussion 
of B. In this same section Sabersky gives a brief survey of some problems in boiling heat transfer. Yuan 
in G discusses film cooling problems from what this author considers to be a rather limited viewpoint. 
In Section I Hottel presents radiation calculations from a practical engineering approach. This section 
will force those who are desirous of any understanding of the actual problems involved to seek further 
information elsewhere. 

Although the volume leaves much to be desired it does nevertheless present a great deal of infor- 


mation which the avid reader can certainly seek out. 
RonaALpD F. PROBSTEIN 


An introduction to mathematical statistics. By H. D. Brunk. Ginn and Company, Boston, 
New York, 1960. xi + 403 pp. $7.00. 


This text is written for students who have studied calculus but have no formal background in 
probability theory or mathematical statistics. The first part of the book is devoted to elements of 
probability theory. A discussion of discrete and continuous random variables follows. The second part 
on statistics is more ambitious. A listing of some or the chapters in this part follows and indicates the 
scope of the discussion: random sampling, the law of large numbers, estimation of parameters, central 
limit theorem, confidence intervals and tests of hypotheses, statistical decision theory, regression, 
experimental design and analysis of variance, distribution free methods. Of course, some of these topics 
are at times treated in an abbreviated manner and without proofs. Nonetheless, the discussion is clear 
and the breadth in coverage startlingly good. Many of the current ideas in statistical practice and 
theory are touched on. There is a very good set of statistical tables at the end of the text. 

M. RosENBLATT 


Exponentially distributed random numbers. By Charles E. Clark and Betty Weber Holz. 
The Johns Hopkins University Press, Baltimore, Maryland, 1960. vii + 249 pp. 
$6.50. 

The book presents a table of 100,000 random numbers from the distribution with probability 
density e~7, 0 < x < o. The generation of these random numbers was based on uniformly distributed 


pseudo-random numbers and the validity of the generating process was tested empirically. The tests 
of randomness and their results, which were satisfactory, are presented in the introduction. 


WALTER F. FREIBERGER 
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Adaptive Control Processes: 


A Guided Tour 


By Richard Bellman 


Dr. Bellman, author and inventor of Dynamic Programming, conducts a guided tout 
through “a few of the methods which work” in the field of adaptive control. His book 
gives a panoramic view of what an ingenious mathematician does when faced with 
the myriad problems of automatic control. The author has minimized detailed rigot 
in the interest of making clear the basic ideas in a broad spectrum of applications. The 
primary purpose is to present a number of techniques which appear promising and 
have already been used with greater or less success. The author takes engineering 
problems which cannot be solved by conventional methods and shows how to get solu 
tions. A RAND Corporation Research Study. 400 pages. $6.50 





Stability in Nonlinear Control Systems 


By A.M. Letov 


Dr. J. P. LaSalle said of this book: “A plain, unsophisticated, painstakingly thorough 
treatise on application of Lyapunov’s direct method. The mathematics required of the 
reader is elementary. The author’s notations . . . possess the advantage (to a person 
wishing to apply results) that the formulas are self-explanatory.” The author, a Nobel 
prizewinner, is a scientist held in the highest esteem by control experts in the U. S. 
He has added to the American translation several additional chapters not included in 


the original. 320 pages. $8.50 





Hydrodynamics 
2nd Edition, Revised and Enlarged 


By Garrett Birkhoff 


\ complete revision of the spirited and stimulating first edition of ten years ago, 
this book brings up to date the results of further activity in this field. The author has 
added a chapter on turbulence, and he has expanded the work on paradoxes and 
modelling. W. M. Elsasser, in the Review of Scientific Instruments, said of the first 
edition: “A book such as this, concentrating as it does on the boundaries of funda 
mental progress, should be indispensable to all those engaged in hydrodynamical re- 
search who are concerned with the type of generalization that so often in the past 
has led to fundamental progress.” 275 pages. $6.50 





Order from your bookstore, or 


PRINCETON UNIVERSITY PRESS 
Princeton, New Jersey 
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